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ABSTRACT 


The Naval Weapons Support Center is planning to implement a bonus system to 
improve the reliability of pyrotechnic devices. The measure of effectiveness that they 
wish to use to determine how to award bonuses is the reliability of pyrotechnic devices. 
The data available to estimate this reliability is based on the current sampling inspection 
plan in which devices are tested in different environments. The models which include 
both dependence and independence assumptions between the outcomes of these tests are 
implemented and estimates of overall reliability along with 95 % lower confidence bound 
are obtained. The 95 % lower confidence bounds are found by bootstrapping. Using 
these estimates, models for making the decision to award bonuses are discussed and 


studied using Monte Carlo simulation . 
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I. INTRODUCTION 


A. BACKGROUND 

Nations spend a lot of money to establish a strong defense network. It is essential 
that nations buy reliable weapons and ammunition from the contractors. To ensure re- 
liability, contracts must include a lot acceptance sampling plan that specifies the minimal 
acceptable quality. 

The Naval Weapons Support Center purchases pyrotechnic devices. Unless other- 
wise specified in the contract, the supplier is responsible to see that his devices meet all 
inspection requirements as specified. The inspection requirements are particular to 
characteristics of each type of device and are specified on the reference drawings and 
supplemental quality assurance provisions of the contract. Any testing that needs to be 
done on these devices is explained in this contract. 

When nations buy weapons and ammunition from the same contractor, they would 
like the quality to improve over time. As contracts are now written, contractors need 
only to satisfy the requirements of the sampling inspection plan for lot acceptance. 
Under such contracts, contractors have no incentive to improve the quality of items they 
provide. For this reason, to improve quality, The Naval Weapons Support Center has 
decided to implement a bonus system. The contractor will be awarded a bonus if the 
result of the sampling inspection exceeds the minimum requirements for lot acceptance. 

The Naval Weapons Support Center will begin to implement a bonus system for 
pyrotechnic devices in FY91. The data available to make the decision whether to award 
a bonus is based on the current sampling inspection plan. This plan is a series of de- 


structive tests in different environments. The purpose of this thesis is to provide the 


Naval Weapons Support Center with guidance for implementing a bonus system for 
pyrotechnic devices. Once implemented, the pyrotechnic bonus system will serve as a 


prototype for bonus systems for other devices. 


B. PYROTECHNIC DEVICE RELIABILITY 

A pyrotechnic device is a chemical and grenade ammunition. There are three cate- 
gories of pyrotechnic devices. The three categories are: Aerial display, Surface display 
and Grenades [Ref. 1: pp. 1-2]. Samples of the Pyrotechnic devices are exposed to 
various environments and then activated. The cniteria for successful activation depends 


on the type of device : 


1. Aerial display (Ground signals, flares, airburst simulators, sub signals, signal kits, 


etc.) 


“Successful activation means that the item will, after simulating user environment, 
successfully .....” 
e launch, 
e have proper separation / signal ignition, 
e reach desired altitude at correct angle, 
e have proper parachute deployment, 
e have proper display color, 
e have proper display time, 
e have no subsequent interference with next item. 
2. Surface (Ground or Water) display (flare, hand-held signals, smoke / illume 
grenades, simulators, etc.) 
“Successful activation means that the item will, after simulating user environment, 


successfully ....” 


e have signal ignition with proper display, 
e have proper display, 
e have proper display time, 


3. Grenades (fragmentation, defensive, white phosphors, etc.) 


“Successful activation means that the item will, after simulating user environment, 
successfully ....” 
e function (high order detonation following delay), 
e have proper dissipation of payload, 
e have completely consumed payload. 
C. PYROTECHNIC DEVICE RELIABILITY PROBLEM 


Samples from any large lot of pyrotechnic devices submitted by a manufacturer must 


activate after exposure to different environments. These are; 
1. Manufacturer Environment, 
2. Temperature and Humidity Environment, 
3. Vibration Environment, 


4. Altitude Environment. 


All items tested are subjected to the manufacturer environment. However items are only 
subjected to one of the three remaining environments : Temperature and Humidity, 
Vibration or Altitude. 

The sampling plan consists of using four distinct samples from a lot that can be as- 
sumed (approximately) statistically independent. The items tested in each sample are 


also assumed to be independent. According to the sampling plan; 


e 20 items are subjected to the Manufacturer Test, 


e 20 items are subjected to both the Temperature and Humidity Test and Man- 


ufacturer Test, 
e 32 items are subjected to both the Vibration Test and Manufacturer Test, 


e 20 items are subjected to both Altitude Test and Manufacturer Test. 


A total of 92 items are tested. 

Acceptance criteria for each test are : 

1. Manufacturer Test : Of the 20 items; if no more than | fails to activate, the lot 
passes. 

2. Joint Temperature and Humidity and Manufacturer Test : Of the 20 items; if no 
more than | fails to activate, the lot passes. 

3. Joint Vibration and Manufacturer Test : Of the 32 items; if no more than 2 fails 
to activate, the lot passes. 

4. Joint Altitude and Manufacturer Test : Of the 20 items; if no more than |! fails to 
activate, the lot passes. 


The number of failures for these tests will be summarized by the vector : 


( FOM, FOTH, FOV ,FO4 ) (1.2) 


where 
¢ FOM represents the number of failures after the manufacturer test, 


e FOTH represents the number of failures after both the temperature and humidity 


and manufacturer tests, 


¢ FOV represents the number of failures after both the vibration and manufacturer 


tests, 


e FOA represents the number of failures after the altitude and manufacturer tests. 


For example, 
(al lees eh). (1.2) 


represents, | failure in manufacturer test, | failure in the joint temperature and humidity 
and manufacturer test, 2 failures in the joint vibration and manufacturer test, and | 
failure in the joint altitude and manufacturer test. It is also the maximum number of 
failures in each test that still leads to lot acceptance. 

The marginal distributions of the number of failures in each of the above tests are 
modeled by binomial distributions. There are 24 possible realizations of the sampling 
inspection; ranging from the best case with no failures to activate, to the worst case with 
1, 1, 2 and 1 devices failing to activate in these tests manufacturer, temperature and 
humidity, vibration and altitude respectively. These cases are tabulated in Table I. 

To award bonuses we need one measure of effectiveness for pyrotechnic devices that 
can be estimated from the available data. Ideally, this measure is the reliability of the 
device. However, because each of the tests are destructive, there is no one natural defi- 
nition of reliability for these devices. To be on the conservative side, we define the reli- 
ability of a device to be the probability that the device will activate after exposure to all 
of the environments. It is not easy to estimate this reliability from the sampling plan 
data. This data is incomplete in the sense that we have limited information about the 
joint probability of activation after exposure to more than one environment. To try to 
compensate for this lack in the data, we will use models for the joint distribution of 
( FOM, FOTH, FOV, FOA ) that specify particular types of dependence between the 


events that a device activates after exposure to different environments. 


Table 1. POSSIBLE CASES 


Using these models and based on sampling plan data, estimates of the overall reli- 
ability along with lower confidence bounds are obtained. These will be used to imple- 
ment the bonus system for pyrotechnic devices. We compute the maximum likelihood 
estimator ( MLE ) of the reliability by maximizing the appropriate likelihood; lower 
confidence bounds ( LCB ) are found by bootstrapping. The MLE’s are computed under 
both independence and dependence assumptions. The estimation procedures assuming 
independence are described in Chapter II. In Chapter III we incorporate dependence 
by fitting a Log Linear Model to our data. The MLE’s from Chapter II and Chapter 
III lead to inappropriate results for this for this problem; thus in Chapter IV we consider 
alternate and very conservative estimates of reliability. Using the estimates of Chapter 
IV, we investigate a sequential scheme for making the decision to award bonuses in 
Chapter V. The results of simulations are presented in Chapter VI. Finally, conclusions 


and recommendations are given in Chapter VII. 


Il. THE ESTIMATION OF THE MAXIMUM LIKELIHOOD ESTIMATOR 


(MLE) OF THE RELIABILITY WITH INDEPENDENCE ASSUMPTION 


A. DEFINITIONS 
We will say that a device survives environment E , if it is still potentially capable of 


activation after exposure to environment E . Let, 
e E, be the device activates after exposure to Manufacturer environment, 


e FE, be the device activates after exposure to Temperature and Humidity environ- 


ment, 
e E, be the device activates after exposure to Vibration environment, 
e £, be the device activates after exposure to Altitude environment. 


We define the reliability of device as below, 


R= P(E 1 1) Ef) £a). (2.1) 


In this formula, R means the probability that a device activates after exposure to four 
environments. We will estimate R for each of the 24 cases which lead to lot acceptance. 

Let QO, = P(E,) be the probability that device activates after exposure to environ- 
ment 1, fori = 1, 2, 3, 4. In the acceptance sampling plan several of the items must 
activate after exposure to a joint manufacturer and another environment. To avoid 
confusion we will denote tests | through 4 as the manufacturer test, the joint temper- 
ature humidity and manufacturer test, the joint vibration and manufacturer test and 


joint altitude and manufacturer test respectively. 


Let 


Re ee) =O, (2.2) 


and let 
R, = P(E) &). (2.3) 


fete x, is the probability that a device survives test i fori = 1, 2, 3, 4. The simplest 
models to assume that £,, £,, £,, £, are independent. If we assume that £, ..., E, 


are independent then 
Ra On; (2.4) 
fori = 2, 3, 4 and the reliability of device is, 
R= QQ) Q3 Qa. (2.5) 


B. THE LIKELIHOOD EQUATION 
eet 
e xX, be the number of devices that activate after test 1, 


¢ n, be the number of items given test 1. 


Then 4X; is binomial with parameters R, and n, for i=1, 2, 3, 4. Under the assumption 


of independence the joint likelihood function of observing X, = x,, ..., X, = 41S 
4 
a x n,— X; 
E(x, , X2,%3,%4| R,,Ro,R3,R) = eel a ha) (2.6) 
= 


i= 


with constraints : 


(= R 21 


OS iG =. 


O<R; 


lA 
Da 


O's ky = Ry 


Our aim is to maximize this likelihood function subject to the constraints that 


(R,, hoo hy, Raye S where 


From the equation (2.6), we see that maximizing L is equivalent to maximizing 


l= ecG: In R;) + (n; — x;) In(1— R;)}, (2.7) 
i=] 


n 


x 


) fori = 1, 2,3, 4 have been dropped ( because they 
do not effect the maximization procedure ) and the natural logarithm of L is taken. 


where the constant multipliers ( 

We first show that / is a concave function. To show that / is a concave function, 
we can show - 7 is a convex function. According to Theorem 3.3.6 [Ref. 2: p. 92], 
by looking at its Hessian matrix, we can learn whether function is convex or not. If its 
Hessian matrix is positive semi-definite at each point S then function / is convex. To 


create the Hessian matrix, we must calculate partial derivatives of the function - /, 


nee in ni — Xj 
oR, = x eae U2) 


~ 7. = — + (2.9) 
aR; j eS 
al ae : 
= as 2.10 
ORR, 0 ix] (2.10) 


Then, the determinant of the Hessian is: 


4 
| H | | |= Fag niet, (2.11) 
R? 1 aaa R? , = 
i=] 


Clearly we can see that forO < R, < 1i1= 1,2,3,4, | H |is always positive. Because 
- / is continuous, this implies that - / is a convex function on S. As a result of this, 
! is a concave function. 

We note that with the constraints on the probabilities R,, ..., R,, there does not 
in general exist a closed form solution to MLE. However, with only 24 realizations of 
X,, ..., X, Of interest, the estimated reliabilities for these 24 cases can be found with 


some rather tedious but straight-forward computations. 


C. COMPUTING THE MAXIMUM LIKELIHOOD ESTIMATORS AND LOWER 
CONFIDENCE BOUNDS 
Because / is concave over the convex set S, if the maximum occurs in the interior 


of S, it is a unique maximum and is given by 
A x; 
R= = (2512) 


fori = 1, 2,3, 4. where 


1] 


> 
lA 
>» > 
lA 
SS 
lA 
D> 
lA 
a > 


for r= 23,4, 


In this case the MLE’s for Q,1 = 1, 2, 3, 4 are 


Oe aie (2.13) 


(2.14) 


for i = 2,3, 4. Finally we can estimate the reliability of the device as 


R = 0, 0,0; Oy. (2.15) 


Table 2 summarizes the cases for which the MLE’s can be found using (2.12) - 
(2.15). When the maximum falls on the boundary of S, there is no explicit expression 
for the MLE of (R,, ..., R,). These are the cases with the exception of ( 1 1 21 ) 
which have one failed item in the manufacturer test. This implies that R, will be less 
than 1.0. To find the MLE, we find the maximum of (R,, .., R,) on each of the 
boundaries, compute / for each of these and let the MLE be the one with the largest 
value of /. It is clear that in most cases several of the boundaries can be eliminated 


from consideration, simplifying computation considerably. 


Table 2. EASY AND HARD CASES 


<eR tee ERCRR SE 


D. EXAMPLES 


1. Easy Case 


In this example our failure vector is, 
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According to this failure vector, we can write our likelihood equation by using equation 


(2.6) to get 


LE = R?(1— Rk PRY (1—R,) Rs (1 — RR, (1- 2) ee 
C= Rk atl, Ue Rea ae 


fori = 2,3,4. From the likelihood above; 


Ry = Sb = 2 = 1.0000 
R, = 3 = B = 0.9500 
R= = = 2 = 09375 
R= = = x = 0.9500. (2.17) 


Results imply that all R's fori = 2, 3, 4 are between 0.0 and R,, which means that 


constraints are met. Now we can estimate Q's using equations (2.13) and (2.14). 


Q, = 1.0000 


Q, = 22300 = 0.9500 


owe oe = 0.9375 (2.18) 


Q,= —2 = 0.9500. 


And finally, reliability of the device can be estimated by using equation (2.15) 


R = (1.0000 ) ( 0.9500 ) ( 0.9375 ) (0.9500) = 0.84609375 . (2.19) 


2. Hard Case 


In this example our failure vector is 


According to this failure vector, we can write our likelihood equation by using equation 


(2.6) as 


oR (l= kyr hel Rs Re (1 — 9 R, ) Ry 1 = R,)' (2.20) 


st = H = 0.9500 

= S = 1.0000 

= = = = 09375 

Ss = 2 = 0.9500. (2.21) 
As you see from the above = > > thus, the MLE does not lie in the interior of S. 


We begin by considering the boundary 


The likelihood equation on this boundary 1s, 


L= RY (1 = Rey (2.22) 
O°s his 
From the likelihood equation, 
R= R= R= R= B = 009565. 


Then value of the likelihood with these estimated R’s 1S 
(is eae = 
Lal 99 eG 99 ) =07 Molin. (2:23) 


For the boundary 


the likelihood equation is, 


L=R(1-R YR (1-8, ) (2.24) 
0 = Ry = 10s eae 
From the likelihood equation, 
R= R= 8 = 3 = 09583 R= 42 = 0.9500, 


and the value of the likelihood with these estimated R’s 1S 


69 69 3 


L = ( 
For the boundary 
R, = R, = ki, SR, 
the likelihood equation is, 
DF Cle) Ra Glvoen,.)- 
(27h is ee = he 


From the likelihood equation, 


R= R= Rh, = 2 = 0.9666 
R= 2 = 0.9375. 


Then value of the likelihood with these estimated R’s iS 


— ( 28 8, 2. y2¢ 30 30, 2 2 _ =! 
lea) aq (ay inl a5) = 0-876 x 10" 


For the boundary 


R, = R23 = Ry, WR < R,, 


the likelihood equation is, 
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69 coe | eye ee 7 
ay) (ay Vag) Coq) = 0-723 x 107. 


(2.25) 


(2.26) 


(2.27) 


ER, Cl Rien 


0 =] RFR, ohe 


co) 
lA 
“a 
lA 


(2.28) 


Clearly L is maximized on the boundary of the constraints (2.28), i.e. on another of the 


boundaries of S thus we can eleminate this case from consideration. For the boundary 


R, = Rk, ,R8 <5 R8,%5 R,, 
the likelihood equation is, 
L = Ry CU 9Ry Rae), Raa Ra) 
OS RS 1, 0S Ry aR, 0 a Ree 
From the likelihood equation, 


Ri = hy == = 09750 Ro 0 oe 


Then value of the likelihood with these estimated R’s iS 


(2.29) 


1 (398; De 17 200 - 


Quick inspection reveals that the remaining boundaries can be eliminated from consid- 


eration. 


After boundary analysis, we can see that equation (2.30) gives us the maximum 


likelihood value. Thus the MLE’s for R,, R,, R;, R, are 


nA 


R, = 0.9750 


R, = 0.9375 
R, = 0.9500. 


Now we can estimate Q,’s with equations (2.13) and (2.14) 


Q, = 0.9750 
_ 0.9750 _ 
Bae 0.9750 foe 
— 0M3TS. 2 (2.31) 
Q3 = 0.9750 ae 
0.9500 _ 
Os = 0.9750 — aes 
And finally reliability of the device can be estimated by using equation (2.15) 
R = (0.9750 ) ( 1.0000 ) (0.9615 ) (0.9744 ) = 0.91346100. (2.32) 


We computed MLE’s with a FORTRAN program MLEA given in Appendix 
A. After 24 replications of this program, we can estimate R, ‘’s for 1 = 1, 2,3, 4 for each 
of the 24 outcomes of the sampling plan that lead to lot acceptance. These are given in 
Table 3, Ris given in Table 4. 

We compute lower confidence bounds for each possible case by bootstrapping. 
The bootstrap can be used to produce approximate confidence intervals in an automatic 
way. There are several ways to set approximate confidence intervals with bootstrapping. 
These are the percentile method, the standard method, the bias-corrected percentile 
method and the nonparametric method [Ref. 3: pp. 67-70}. Let 6 be an unknown pa- 


rameter with estimator 6. To bootstrap, samples are generated using @ in place of the 
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unknown parameter 8. For each of these samples an estimator 0° is computed. We 
define G (s) to be the parametric bootstrap cumulative distribution function of 8", Le 
G(s) is the empirical distribution of the 6’s. All methods mentioned above use 
percentiles of G to define the confidence interval. They differ in which percentiles are 
used. The percentile method was used in our calculations. 

lf we use the notation @ (a) for the level « endpoint of one sided lower confi- 


dence interval for @ , then 
PLO(a) = oaie=eor (2:33) 
An estimate of 8 (a) from the bootstrap cumulative distribution is given by 
A Ny) 
O(a) =G (a). (2.34) 


To get LCB’s for R, for each of the 24 realizations of the sampling plan , we 
generate bootstrap samples of random failure vectors, in which failures come from in- 
dependent binomial distribution with parameters (17,, R,) fori = 1, 2,3, 4. This is done 
using the FORTRAN program RANVEC [Ref. 4] in Appendix B. We generate 5000 
failure vectors for each case. And then we estimate R ’s from each of the 5000 failure 
vectors by the means of the program MLEA. The next step is to compute the order 
statistics of R ‘s from R, to Rae We get the parametric bootstrap cumulative distrib- 
ution function ( i.e. the emprical distribution of R, eee Re ) with this computation. 
This is done using the FORTRAN program SORT in Appendix C. 

Finally we can get the 95 % lower confidence bound using equation (2.34) from 
this routine. Reliabilities and 95 % lower confidence bounds are listed in Table 4. Re- 
sults in Table 4 are given in descending order and R’s and LCB’s are not ordered as we 


expected them to be. For example, the failure vector that has the maximum number of 
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failure in each test is in the middle of the table with respect to R. Ithas also a bigger 
95 % LCB than the failure vector that has a total of 3 failed items in each of the joint 
tests. These results are counter-intuitive because we expect that more failures indicate 
a lower overall reliability. This is reasonable because it 1s likely that a device that is 
poorly constructed is more likely to fail in any of the four environments. From the re- 
sults in Table 4, it is clear that an attempt to model dependence must be made in order 


to get believable estimates of R. 
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Table 3. TEST PROBABILITIES ( A ) 
[<ASE | ssa a 
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Table 4. RELIABILITIES AND 95 % LOWER CONFIDENCE BOUNDS 
(A) 


T0000) 
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II. LOG LINEAR MODEL WITH DEPENDENCE ASSUMPTION 


A. BACKGROUND 

Modeling the outcomes of the various environmental tests as independent is clearly 
inappropriate. However, the nature of the acceptance sampling plan makes it impossible 
to estimate the reliability of the device ( i.e. the probability that it would activate after 
exposure to all four environments ) without some assumptions about the dependence 
between outcomes of various tests. Thus, our second approach is to model the 
pyrotechnic device reliability with a log linear model. 

The results of our test series create a (2x 2x2x2) contingency table. Let p 


represent passing a test and f represent failing a test. 


Table 5. CONTINGENCY TABLE STRUCTURE FOR TEST SERIES 


Passed Failed 
Manufacturer Test Manufacturer Test 


Passed Failed Passed Failed 
Tem.&Hum. | Tem.&Hum. | Tem.&Hum. | Tem.& Hum. 
Test Test Test Test 
Passed 
Altitude Moppp Mo fpp Myppp Mrpp 
Test 
Failed 
Altitude Moppf Mo spp Mypps Moor 
Test 
Passed 
Test 
Failed 
Altitude Mop Mor M psy Mpy 
Test 
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Passed 
Vibration 
Test 


Failed 
Vibration 
Test 


The results of our tests can be thought of as censored data from a hypothetical 
(2X 2X2 X 2) contingency table (See Table 1) . The frequency in each cell of this table 


is 
Mjyr , (reais k, lve Coes 


the number of devices out of 92 which would have result 1 in environment | (manufac- 
turer), result j in environment 2 (temperature and humidity alone), result k in environ- 
ment 3 (vibration alone), result / in environment 4 (altitude alone). 

This is a hypothetical table, because if a device was exposed to all four environments 
and then failed to activate, there would be no way to discern which combination of the 
four environments caused failure. The data from the acceptance sampling plan can be 
thought of as censored data from such a (2 X 2 X 2 X 2) contingency table. As an 
example, a device that is given just the manufacturer test belongs in one of the cells 
(p,j,k,/) where (j,k, /) e {p, f}3, because it is not clear what would have happened 
to it had it been exposed to the other three environments. 

Using a log linear model under independence, the expected value of each cell fre- 


quency 1s 
1 2 3 4 
E[ Miu J ae ee ee Ae ey (3.1) 
where i, j, k, / = p,fand 


—— ry (3.2) 


fori = 1, 2,3, 4. Here the parameter yu represents the overall effect and the parameters 
4‘ represent the effect of passing environment i. It is simple to show that the log linear 


model (3.1) is equivalent to the independence of the environments [Ref. 5: pp. 25-46]. 
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On a logarithmic scale, the independence relation is equivalent to the additive re- 


lationship. As an example, 


log M 


pppp = H+ dp zis IE ale A, * ie (E58) 


Log linear models, which take into account dependence include extra interaction terms. 
Because of the extreme amount of censoring, we will only consider models with two way 


interaction terms of the type. 


1 2 3 4 12 13 14 23 24 % 


where i,j,k, / = p, f, 


A= -% (3.5) 
fori = I, 2, 3, 4 and 
ih = tym = a= - ah 06 
for i,j, = 1, 2, 3,4. As an example; 
PIL Moppp J = gl Apt Ay ay ty Ape et tee tye tp (3.7) 


The censored data does not tell us much about partial association between pairs of 
tests. For this reason, we assume that partial associations between any pair of tests are 
all the same. Let @ be the partial association between tests. Then we can reformulate 


equation (3.7) as follows: 


i 2 3 4 
ENG 552 — gta, +4, tata, +60 (3.8) 
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The MLE’s for the expected number of devices in each cell ( or equivalently the 
MLE’s for the parameters ) can not be found explicitly. The EM Algorithm ( Expecta- 
tion - Maximization ) will be used to approximate the MLE’s of the expected number 


of devices in each cell. 


B. EXPECTATION MAXIMIZATION ( EM ) ALGORITHM 

A general method of maximum likelihood estimation from incomplete data is the 
EM Algorithm. The expectation maximization algorithm 1s an iterative procedure where 
each stage consists of: 

® an expectation step ( E ) followed by 

® a maximization step ( M ). 

This algorithm is generally used to compute maximum likelihood estimators in in- 


complete data problems. In the application of the EM algorithm we: 
e replace missing values by their estimated expected values given the incomplete data 
¢ estimate parameters 
e reestimate the missing values assuming the new parameter estimates are correct 
e reestimate parameters 


and so forth, iterating until convergence [Ref. 6: pp. 127-141]. This iterative algorithm 
works as follows in pyrotechnic device problem. 
1. Initialization 
The algorithm requires initial guesses for the parameters of the log linear model 
cell frequency in contingency table. It uses initial guesses and calculates initial cell 
probabilities. We will get our initial guesses by first fitting a log linear model under in- 


dependence. 


27 


2. Iterations 
= Expectation Step 
EM Algorithm estimates expected cell frequencies given the data by using the 
most current estimates of the parameters. It compares these estimated expected cell 
frequencies with the previous estimates. If the differences between these two estimates 
are small enough, then the EM Algorithm has converged. We then accept the final es- 
timates of cell frequencies as MLE, so we can easily calculate cell probabilities. 
s Maximization Step 
If the algorithm has not converged, then it starts to estimate new parameters for 
the log linear model using expected cell frequencies from the E step as if there were ac- 
tual data available. Estimation is done by maximizing the likelihood using an iterative 
Newton-Raphson method. The estimated parameters from this step are than used in the 


next E step of the EM algorithm. 


C. CALCULATIONS 

We apply the EM algorithm for each realization of the failure vector ( FOM, 
FOTH, FOV, FOA ). During the maximization step of the EM Algorithm, we use the 
Newton-Raphson procedure which is described by SAS, [Ref. 7: pp. 190-212]. We first 


describe the M step using as an example the following failure vector. 
(. Dae 2ao) 


s Maximization Step 
We start with initial guesses (IG) for M,,,s. These are the conditional expected 


values provided by the previous E step. For example, initial guesses for the failure vector 


above, are shown below. 
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Pppp 
Mppps 
Moptp 83.08 
0.1 
Mopiy 3.36 
Mpyppp 0.1 
Mpypr 2.21 
M 0.1 
7 0.1 
Moy 0.1 
IG = [ore ea (3.9) 
Mrppp i ‘i 
M : 
Sppef 0.1 
M ppp 0.1 
My 0.1 
Mpp 0.1 
0.1 
Mor 0.1 
Mrpp 
Mypy 


Mppf 0.9030 
92 0.0011 


PPsp = 0.0365 : (3.10) 


~~> 
lI 
x 


0.0011 


We will use P to get MLE’s for the parameters of the log linear model. Our log linear 


model which includes two way interactions can be written as 
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log P 


PPpp 
] 1 6 
18 vprl : 1 Stee 
log Popyp = eee L (3.11) 
. 1 -1 -1 -1 —1 6 ; 
log Pry 


where y7 = (yu, J!, 42, 43, 44) is the parameter vector and p is normalizing constant 
required by the restriction that the probabilities sum to 1 and where P,,,, is the proba- 


bility corresponding to the (ijk /)* cell in the hypothetical contingency table. We can 


derive from equation (3.11) 


log Prppp 

log P de se Sie 
is Lt jie 10 

log Prop = ee in) |e (3.112) 
. —-1 -—-1 -1 -1 6 , 

log Pry 


where B7 = (J!, 42, 43, A$) is the parameter vector and 6 is normalizing constant re- 
quired by the restriction that the probabilities sum to 1. Then, to use the SAS proce- 


dure, we rewrite log P,,,,, 


log P 


pppfr °° 


. , log Pay as the 15 logits, F, = (log P,,,,/ log Pry), 


., Fs = (log Pp! log Pyy) so that 


F log Poppp 
F, 0140 0 =I digs 

F(P)=] - = Pre a log Prorp (3.13) 
in ONO! 0 0 =e 
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where P = (P,,,,, Pon +--+ Pry). Using equations above the following result is ob- 


tained 


a4 
Non 
NO ho 
t 
No 
=) 


2) ol ee eet ee (Ee (3.14) 
es eae 


The design matrix (X) of this problem, from the equation (3.14) is as follows: 


(3.15) 


be 
] 

SOO CO COCO CO ON NNN NWN WN 

SOON NNN DCO ON WN WN VY 

SCNnNndNIGDONN CONN CON VY 

NON ON ON OCON ON OCN OW 
A 


In the application of Newton-Raphson Method, we use the variance and covariance 


matrix S of F( P). It’s inverse is given by 


pepe ee — Poop * Popp ee x Pigp 
= _ p oe aR 
s7(P) a "pepe * Proppf —Ppppf P ppt hae ug Site (3.16) 
a 2 
z Popp “ Pip © gd Pip Pri a Pint 
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The first estimates of parameters for the Newton-Raphson procedure are calculated 


as follows: 
b, = [Xess CP)X PL X° Sees a (3.17) 


where S-' and F are estimated using the proportions in the vector P. We estimate the 


reduced logit response functions using the equation 


= by, (3.18) 
From F, we can compute the updated estimates of Po, Po, ---5 f),0eee 
Il = (T1,(1), 1,(1), ..., T,.(1)) be the vector, which contains these estimatesion 


the cell probabilities. The value of the log likelihood evaluated at I is 


16 
LHE = ) x; log 11(/) (3.19) 
i=1 
where x, is the number of items in cell i with respect to probability in T1. 
We estimate parameters iteratively until the difference between last estimate and 
previous is small enough. At each iteration we update the inverse of variance and 
covariance matrix with probabilities of II from the previous iteration and we do fol- 


lowing matrix computations. 


C=205- “eller 


, 3.20 
G = X'[920 x (P-—N)] Cae 


Let 6, be the next estimate of parameters £ in the i th iteration. Then 4, as follows: 


b, = b_, - 6C'G (3.21) 


o2 


where 6 < | is a constant supplied by the user. 

We get first cell expectations from initialization. The failure vector, which is defined 
at page 5, has initial guesses for cells which are shown in equation (3.9). After the first 
application of Newton -Raphson Algorithm, we have following initial cell expectations 


(CE ). 


82.4278259 
0.7278117 
3.1047230 
0.1208856 
2.0790358 
0.0809494 
0.3453168 
0.0592893 
Oe Neco Cz) 
0.0789160 
0.3366428 
0.0578000 
0.2254283 
0.0387050 
0.1651089 
0.1250073 


= Expectation Step 
The conditional expected frequency for each cell expectations are calculated using 
the estimated parameters of the log linear model from the previous M step. Some of the 


cell expectations formulated below 


ELM, 


ppp | (FOM, FOTH, FOV, FOA) | = (20 — FOM)x P(pppp|p...) 


+ (20 — FOTH)x Pp ppeleP-») 553 
+ (32 — FOV)x P(pppp\p.-p.-) 
+ (20 — FOA)x P(pppp|\p.-p) 
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E UMpyror | (FOM, FOTH,FOV, FOA)] = (20 — FOM) x (pi pi | pas 
+ (FOURS py psf .-) (3.24) 
(32 ON a cpt pif |i peste) 
+ (FO ADPAa ie \p.-J) 
where ” . ” represents either a pass or fail. For example, 


P(pppp|p...)is the probability that device passes all environments given that it 
passed manufacturer test. 

P(pppp\|pp.-)is the probability that a device passes all environments given that 
it passed manufacturer test and temperature-humidity environment. 

The conditional probabilities are computed from the estimated cell probabilities 


from the previous M step. As an example, we have following failure vector, 


em. Sau 


Then EL M,,,|(1, 1, 2, 0) ] is calculated as follows: 


PPPP 


| 82.4278259 
P(pppplp..-) = 88.9458375. = 0.9267193 


82.4278259 
Pipppplpp..) = "86.3812462. 0.9542328 


82.4278259 
P(pppp|p.p.)— 85.3156228 0.9661516 


82.4278259 
P(pppplp ep? = 87.9569015. = 0.9371388, 


where the numbers come from (..22) the estimated cell expected values after from the 


previous M step. Thus, 
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EL Mpppp | (1, 1, 2, 0) J = 83.4654163 


The remaining 16 conditional expectations are calculated similarly. Then we com- 
pare these expectations with previous expectations. If the difference between compared 
expectations are small enough then EM Algorithm is assumed to have converged. 

When EM Algorithm converges, we can estimate reliability of device ( R ) after 
exposure all four environments as 

A E [ M LM pppp 4d ] 


a ee (3.25) 


For example, for the failure vector above, after 27 iterations, we obtain 


Bll Myopp | = 84.412714, 
and 
® _ 84412714 _ 
Fe eee WUE 


D. INITIAL GUESS PROBLEM 

The EM Algorithm uses an initial guess vector as in equation (3.9). We use the 
independence assumption as described in equations (3.1), (3.2) and (3.3). The following 
procedures are used during the application of EM Algorithm to the log linear model with 
the independence assumption. 

We have still an initial guesses problem for parameters in the log linear model with 
the independence assumption. For this reason, we do some calculations to find initial 
guesses for the parameters. Moreover, one can choose random initial guesses for the 


parameters too. 
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Assume that we have the following failure vector, 


In this procedure we divide the number of devices which pass the manufacturer test 
by the number of cells which include items that have passed the manufacturer test. In 


this example, we can calculate M,,,, as follows: 


= 2.375 items for manufacturer test 


0.125 item for vibration test 


= 
I 


19.0 
8.0 
[ 20.0 J] = 2.500 dj 
oo items for temperature and humidity test 
1.0 
0 
0.0 


0.000 item for altitude test 


® 
al 
ii 
Lo 
I 


finally the expected number of device for this cell is, 


My = 5.000. 
After computing M,,, for each cell using the same procedure above, we can group these 
cells to estimate Mo, M,oo, Myw, Myo as follows: 
© Mau = Myo 
© Mn = Mp0 
Man = Mn 
© Ma = Mao 


where k, / = p,f. 
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In this example, 


Table 6. GRUOPED DATA FOR INITAL GUESS CALCULATION 


TERMS ESTIMATED VALUES 
1, 00 


By combining the parameters » and /!, we have four equations and four unknown 


parameters. We can easily solve for the parameters from the initial guesses. 


Table 7. ESTIMATED PARAMETERS FOR INITIAL GUESS 


Finally, we can split the sum of the y» and 4! into two parts. One possibility is to 
divide the sum by two assigning half to » and half to 2!. Results using a random initial 
guess and the above procedure for an initial guess are very close. Calculation of initial 
guess parameters are done by the FORTRAN program INITIAL in Appendix D and 


the FORTRAN program PARAM in Appendix E. 
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E. RESULTS 

Reliabilities of the pyrotechnic device are calculated by a FORTRAN program 
LLMDEP in Appendix F. They are given in Table 8. Results from the log linear model 
design with dependence assumption are similar to the first model. We can easily see 
from the previous table that the ordering of the estimates of R is counter-intuitive. We 
expect (1 1 2 1 ) to yield the smallest R rather than an R = 0.934 which is larger than 
R for about half of the 24 realizations of the failure vector. 

It is clear from these results that this log linear model is inappropriate for modeling 
the outcomes of the sampling inspection plan. More realistic models would include 
three-way and four-way interaction terms. However, due to the extreme censoring in the 
data, we can not estimate R for these models. Moreover, reliabilities of several cases for 
even this model were not calculated because of computational limitations. Taking the 
account of dependence with log linear models is clearly not a reasonable approach for 


estimating reliabilities from this data. 
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Table 8. RELIABILITIES WITH LOGLINEAR MODEL 


[fanueevecton | A 


oy 


IV. WORST CASE SCENARIO 


A. ASSOCIATION ANALYSIS 
In the inspection sampling plan, failures are not assigned to an individual failure 
mechanism except for the manufacturer test. For example, from the following failure 


vector. 


We know that there is at least one failure due to the manufacturing test. There is one 
failure from the joint temperature and humidity and manufacturer test, but we do not 
know which failure mechanism generated this failure. Items can fail due to a manufac- 
turing related failure mechanism or other failure mechanism or both. The worst case 1s 
to assume that the cause of failure is due to all of the failure mechanisms that it was 
exposed to. 

Let MT represents manufacturer test, TAHT represents temperature and humidity 
test, VT represents vibration test, AT represents altitude test. 

Because all 92 items are exposed to manufacturer test, the worst case of the failure 


analysis of the failure vector above is shown by the following table. 


Table 9. FAILURE ANALYSIS 


# OF SUCCESFUL 
ITEMS 


This ” worst case ” scenario should give us lower bounds for the estimates R of R. 
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If we calculate the reliability of series system assuming failure mechanisms in test 
are independent when in fact they are associated but not independent, then we underes- 
timate system reliability. It is reasonable to expect that there 1s positive dependence 
between four tests. If one item fails in manufacturer test, it is more likely to fail in other 
tests. One way of modelling positive dependance between tests 1s to assume that tests 
are positively quadrant dependent. 

Given random variables 7,, ..., 7,. They are said to be positively quadrant de- 


pendent (PQD) [Ref. 8: p. 33], if 


BCR See Sl Gre (4.1) 
=1 
for all (4, 4, ..., ¢,) € R«. An equivalent formulation of positive quadrant depend- 
Cieesis 1, 15, ..-, are PQD: iff 
: n 
(Disp Saree = tie ||| /2en ea) (4.2) 
i=) 


for all (4, 4, ..., t,) € R* The proof that (4.1) and (4.2) are equivalent as given in 
[Ref. 8: pp. 32-33]. 
We may take account positive dependence between tests with the PQD assumption. 


Let 


1, if an item passes environment i 
; (4.3) 


0, if an item fails environment i 
fori = 1, 2,3, 4. Assume that 7,, 7,, 7T;, T, are positively quadrant independent, i.e. 
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4 
P(T, 2 5 hee <5, 7, Seem ltt =) (4.4) 
i=} 


or equivalently that 
4 
LUPE ne ere nS | PCRS ah (4.5) 
i=1 


Let R be the probability that an item activates after exposure to all 4 environments. 


Then 
R= P(t, =1, GG =), = 1 (4.6) 
Using equations (4.4) and (4.5) 
Roe P(T, = 1) P(T, = 1) PUT, = 1) PCT] = 1). (4.7) 
Using the notation from the previous section 


R= Q00;9%, (78) 


where Q, is the probability that an item passes environment i. With the censored data, 
we can’t estimate Q,, Q,, Q,andQ, without building a more structured model for 
y re ee acy aval 18: One alternative is get a lower bound for estimates of 


Q,, 2,, Q;, Q,. This is worst case scenario . 


B. CALCULATIONS WITH EXAMPLE 
Let 
» R be the MLE for the reliability of the device. 


s QO, be the MLE for the probability that an item passes environment i. 
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= R bea lower bound of R for reliability of the device. 


= Q, be a lower bound for the probability that an item passes environment i. 


Then according to equation (4.8) 


Further we will construct estimates Q, of Q, such that 
QO; < Q; 


and then define 


»> 
lI 
cate 
tO> 


Thus 


>> 
IA 
os) 
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(4.9) 


(4.10) 


(4.11) 


(4.12) 


From Chapter II, the likelihood of observing X, = x, X, = x, X; = 5, X, = x, 18 


4 


Nn; = 
L (x; .%, Xa %4) Ri oe a | |( ) Ra SR 


x; 
i=] 


- (") iGreen sme 


xy 
Ry 
x P(T, =1,T =1)? 
x2 
(P(t = — 1) eee (4.13) 


x 


13 . 
P(T, =1,7,=1)° 
X3 


ye P(t) — let, de ee 


Mg 
x ( )pa- uno" 


X%4 


¥ 1 P= 1 y= 1) 


where x, is the number of devices out of m, that activate after testi. If we know why the 
device failed for the tests i = 2, 3, 4 which include manufacturer test along with exposure 


to environment 1 then our likelihood could be written as 


ny 
Pa Pr =o" 


xy 


Ny 
x 
X25 X21, X22, 23 


PR = 1, 4 =1)7 PT, = 1, GB =0) 


L(x, %2, Xj, X41 Ry, Ro, Ry, Ry) = ( 


Xx 


Xx 


P(T, =0, T) = 1Y” P(T, =0, T) = 0) 


ny 
(4.14) 
X31 X319 X32. 433 


P(T, =1, T) = 1)" P(T, = 1, T; = 0) 


x 


XxX 


P(T, =0, Ty = 1) P(T, =0, T, = 0) 


( ‘ 
X41 Xai» X42, X43 


P(T, =1,T=1)* P(T, =1, % =0)" 


XxX 


x 


XxX 


P(T, =0, Ty = 1%" P(T, =0, T, = 0)" 


Xx 


where x, 1s the number given test i which failed due to environment i but passed manu- 
facturing, x, is the number given test i that passed environment i but failed manufac- 
turing and x, 1s the number given test 1 that failed both manufacturing and environment 
ifori = 2,3, 4. Note that n, — x, = x, + x2 + x, fori = 2, 3,4. 
From this likelihood, the MLE’s of Q, = P(T,=1) are given by 
= Xp + (XQ +1) + (3 +431) + Oa + Xai) 


LS a 9 ne (4.15) 
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nS ae (4.16) 
0,=- 3 z 7 (4.17) 
Oe ae (4.18) 
However since we do not know the x,’s, we see that 
é = eal aad nL (4.19) 
and 
2 He. (4.20) 
fori = 2,3,4. Let 
5, —_ = =f ae X3 + X% (4.21) 
and 
0= = (4.22) 
fori = 2,3,4. Then O's are lower bounds for the true MLE’s or and 
R i 0, QO, 0; Or (4.23) 


It is clear that R is a lower bound for the MLE R. 
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In Worst Case Scenario, we assume that the item, which fails in test 1, fails because 
of both manufacturing related failure mechanism and other failure mechanism which is 


induced by the i th test environment. With this assumption, the reliability of of device 


is given as follows 


Ro 0,0, 0,0, (4.24) 
where 
A = x; ate x2 oF X3 =e X4 
?; —- 9? 
A x4 
Q, = Sp, 
20 
i x; (4.25) 
2; a 2%) 
A X4 
o 20 
Here is an example, in this example our failure vector is, 
C20. lec al) 
0, Begs har ee 0.945651 
92 
0, = — = 0.9500000 
20 
a 30 : (4.26) 
Q; = oo ae 0.9375000 
ee = 90.9500000 
20 


And finally reliability of the device can be estimated by using equation (4.24) 


R = (0.9456521 ) ( 0.9500000 ) ( 0.9375000 ) ( 0.9500000 ) = 0.8093069 (4.27) 
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C. RESULTS 

We compute MLE’s with a FORTRAN program MLEB in Appendix G. These 
are given Tables 10 and 11. 

To get an approximate lower confidence bounds for R using the worst case data, 
we bootstrap using the procedure described in the previous chapter. The FORTRAN 
program RANVEC in Appendix B is used to generate the 5000 bootstrap samples for 
each case. We then estimate R ’s for 5000 failure vectors by the means of MLEB for 
given case. The next step is to compute the order statistics of R’s from R, to Raw This 
is done by FORTRAN program SORT in Appendix C. Finally, we obtain the 95 % 
lower confidence bound from this routine. Reliabilities and 95 % lower confidence 
bounds are tabulated in the following pages. Results of reliabilities and 95 % lower 
confidence bounds are given in descending order. 

These estimates and LCB’s for R are decreasing with the total number of the failures 
out of 92 items tested. This, at least, is consistent with how we believe pyrotechnic de- 
vices behave. It should be noted that these estimates are in fact conservative lower 
bounds for the time MLE’s under the very weak assumption of PQD. How conservative 
these estimates are cannot be determined without more extensive data that allows us to 


estimate the degree of dependence between tests. 
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Table 10. TEST PROBABILITIES (B ) 
[eS A 
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Table 11. RELAIBILITIES AND 95 % LOWER CONFIDENCE 
BOUNDS (B ) 


[tee | 4d vecron_ | %%1C% 
T0000) 


Ps 


tr 
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V. BONUS SYSTEM APPROACHES 


A. BACKGROUND 

Quality is described as ” especially high degree of goodness or worth ” [Ref. 9 
: p. 685]. In industry, a quality product is one that fulfills customer expectations. 
There are two general aspects of quality: 

® Quality of Design 

® Quality of Conformance. 

All goods and services are produced in various grades or levels of quality. These 
variations in grades or levels of quality are intentional; therefore the appropriate tech- 
nical term is quality of design. The quality of conformance is how well the product 
conforms to the specifications and tolerances required by the design. Quality of con- 
formance is influenced by a number of the following factors: 

® the choice of manufacturing process 

® the training and supervision of the workforce 

® type of quality assurance system ( process controls, tests, inspections, etc. ) used 

= the extent to which these quality-assurance procedures are followed 

= the motivation of workforce to achieve quality. 

Quality Control is the engineering and management activity by which we measure 
the quality characteristics of a product, comparing them with specifications or require- 
ments and taking appropriate remedial action whenever there is a difference between the 
actual performance and the standard [Ref 10: pp. 1-3]. 

As contracts are now written, contractors need only to satisfy the requirments of the 
sampling inspection plan for lot acceptance. Contractors have no incentive to improve 


the quality of the items they provide, although they are in a position to do so. As 


eal 


mentioned before, the quality of conformance is influenced by the motivation. Therefore 
we can motivate manufacturers by giving a bonus for improved quality. To improve 
quality, The Naval Weapons Support Center has decided to implement a bonus system 


for pyrotechnic devices. 


B. BONUS PLANS 

In this chapter, we design a bonus system to improve the quality of pyrotechnic 
devices. A good bonus system encourages reliability growth, because firms try to reach 
a high quality to get a bonus. An effective bonus system must detect small differences 
among the offered lots. In the previous chapters, we estimated the reliability of a 
pyrotechnic device in several ways. We assumed independence in the first model, and 
we assumed dependence between tests in the last two models. The estimated reliabilities 
in the first two models are close to each other, but they exhibit different structures in 
order. One can easily see that the failure vector of the worst case ( 1 1 2 1 ) has bigger 
reliability value than ten of the possible cases, from Table 5 and Table 15. In addition, 
the order of 95 % LCB’s of cases does not match to the order of MLE’s in independent 
models. But in the worst case scenario, Which assumes dependence between tests, we 
get the same order for both 95 % LCB’s and MLE’s. Thus, we will use the worst case 
scenario model in further calculations for the bonus system. With the assumptions 
above, we can apply three sampling plans for giving bonuses to the manufacturers. 
J hese ane: 

# Single Sampling Bonus System 

=* Double Sampling Bonus System 

# Multi-Sampling Bonus System. 

Manufacturers have to meet pyrotechnic device acceptance criteria first, before 


having a chance to get a bonus. 


D2 


1. Single Sampling Bonus System 
In this bonus system, we will decide two things at the end of inspection. First 
we will decide whether the offered lot is acceptable or not, and then for lots which are 
accepted we will decide whether to give a bonus to the manufacturer. 


The Single sampling bonus system 1s designated by three numbers. These are, 


Mio CCOLD. LOD 


where 
= 7, is the sample size for inspection ( here n, = 92 ). 
= LCBFB is the cut off value applied to the lower confidence bound for 
awarding a bonus. 
= LCB means estimated lower confidence bound after inspection. (If LCB is 
greater than or equal to LCBFB, then a bonus is awarded.) 
The following algorithm shows us how single sampling bonus system works for 


pyrotechnic devices. 


STEP # 1 : DETERMINE CUT OFF VALUE FOR AWARDING THE BONUS. 
STEP # 2 : TAKE A SAMPLE SIZE OF 92 FROM OFFERED LOT. 

STEP # 3  : APPLY MANUFACTURER AND THREE ENVIRONMENT TESTS. 
STEP # 4 : COMPARE RESULTS OF TESTS WITH ACCEPTANCE CRITERIA. 
STEP #5 : ESTIMATE ITS LOWER CONFIDENCE BOUND ( LCB ), 


IF THEY MET ACCEPTANCE CRITERIA. 
STEP # 6 : COMPARE LCB OF LOT WITH LCBFB 
STEP # 7 : IF THE LCB IS GREATER THAN OR EQUAL TO LCBFB, 


GIVE BONUS TO THE FIRM. 
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According to the algorithm described above, the following are possible events 
for the firms. 

« Firm may not satisfy our acceptance criteria. This means that firm gets a 
failure vector worse than (1121 ). 

« Firm satisfies acceptance criteria, but its LCB may be less than LCBFB. This 
means that firm does not get bonus, but the lot is accepted. 

« Firm satisfies acceptance criteria, and its LCB may be greater than or equal 
to LCBFB. This means that firm gets the bonus. 

2. Double Sampling Bonus System 

A double sampling plan has an advantage over a single sampling plan. Because 
a double sampling plan involves a larger sample size, it reduces the chance that a man- 
ufacturer who deserves a bonus will not get one. Double Sampling Bonus System per- 
mits the taking of two samples on which to make a decision [Ref. 11: pp. 184-185]. 

In this system, we have two inspection stages. If the firm does not get a bonus 
after the first inspection, then the firm is given a second chance with a second inspection. 


A double sampling bonus system is designated by five numbers. 


ny, n,, LCBFB, LCB1, LCB2 


where 

= 1, 1s the sample size for first inspection, ( here n, = 92 ) 

# ,1S the sample size for second inspection, ( here m, = 92 ) 

« LCBFB is the cut off value applied to lower confidence bounds for awarding 
bonuses, 

= LCB1 is the estimated lower confidence bound after first inspection, 

» LCB2 1s the updated estimate of lower confidence bound after second in- 


spection. 


In the double sampling bonus system, if the firm does not get a bonus after the 
first inspection but does meet the acceptance criteria for the first sample, then a second 
sample is taken. We already have LCB’s of possible cases after the first inspection in 
Table 11. After the second inspection, we calculate the LCB using an aggregated failure 
vector which includes failures from both samples. When we compute LCB’s using ag- 
gregated failure vector after the second inspection, we have a total sample size of 184; 
there are 135 different failure vectors for which the lot meets the acceptance criteria for 
both samples. After tabulating these possible 135 cases and the estimates of R, we used 
the bootstrap procedure to find LCB’s for each case. We created 5000 random failure 
vectors for each of them by using case success probabilities in tests. After this, we esti- 
mated MLE’s of 5000 failure vectors for each possible case. Finally, we estimated 95 
% LCB of each case. These calculations were done by using programs RANVEC, 
MLEB, SORT in appendix B, G and C respectively. The results are tabulated at Ap- 
pendix I. The following algorithm shows us how double sampling bonus system works 


for pyrotechnic device. 


STEP # 1 : DETERMINE LCBFB. 
STEP # 2 : TAKE A SAMPLE SIZE OF 92 FROM OFFERED LOT. 
STEP # 3 : APPLY MANUFACTURER AND THREE ENVIRONMENT TESTS. 
STEP # 4 : COMPARE RESULTS OF TESTS WITH ACCEPTANCE CRITERIA. 
STEP # 5 : ESTIMATE ITS LOWER CONFIDENCE BOUND ( LCB1 ), 
IF THEY MET ACCEPTANCE CRITERIA. 
STEP # 6 : COMPARE LCB1 OF LOT WITH LCBFB 
STEP # 7 : IF THE LCB1 IS GREATER THAN OR EQUAL TO LCBFB, 
GIVE BONUS TO THE FIRM. 


STEP # 8: IF THE LCB1 IS LESS THAN LCBFB, 
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GIVE A SECOND CHANCE TO THE FIRM FOR BONUS. 
STEP # 9 : TAKE A NEW SAMPLE SIZE OF 92 FROM OFFERED LOT. 
STEP # 10 : APPLY MANUFACTURER AND THREE ENVIRONMENT TESTS. 
STEP # 11 : AFTER GETTING THE NEW FAILURE VECTOR, 
ADD THIS ONE TO THE FIRST FAILURE VECTOR. 
STEP # 12 : ESTIMATE ITS LOWER CONFIDENCE BOUND ( LCB2 ), 
WITH AGGREGATED FAILURE VECTOR. 
STEP # 13 : COMPARE LCB2 OF LOT WITH LCBFB 
STEP # 14 : IF THE LCB2 IS GREATER THAN OR EQUAL TO LCBFB, 


GIVE BONUS TO THE FIRM. 


According to the algorithm described above, the following are possible events 
for the firms. 
= The firm may not satisfy our acceptance criteria after first inspection. This 
means that firm gets a failure vector worse than (1121). The lot is not accepted and 
second sample is not taken. 
= The firm satisfies acceptance criteria, and its LCBI may be greater than or 
equal to LCBFB. This means that it gets bonus after first inspection, and a second 
sample is not needed. 
= The firm satisfies acceptance criteria, but its LCB1 may be less than LCBFB. 
This means that it does not get the bonus after first inspection, but still has a chance to 
get a bonus if it submits a second sample. 
Of those which submit a second sample 
= LCB2 may be less than LCBFB. It means that firm does not get the bonus, 


but the lot is accepted. 
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= The aggregated failure vector meets acceptance criteria. But LCB2 may be 

grater than or equal to LCBFB. It means that firm gets bonus. 
3. Multi-Sampling Bonus System 

In this system approach, we may take a sample of size 92 from offered lot three 
times, four times or more. But if we use the three sampling bonus system, after the third 
inspection, We may get(4 X4X7x4 = 448 ) 448 different failure vectors. Multi- 
Sampling Bonus System is going to be more computationally intensive with respect to 
double sampling bonus system. This can be done later using with the same reasoning 


in double sampling bonus system. 


C. EXAMPLES 
1. Single Sampling Bonus System 
Let us assume that we decided that lower confidence bound for bonus will be 
0.8500. We are going to use this LCBFB in these examples. 
= Case # 1 


At the end of inspection, firm has following failure vector 
(42-2020) hs. ji 


The firm did not meet the acceptance criteria. We immediately reject the lot and no 


bonus 1s given. 
= Case # 2 


At the end of inspection, firm has following failure vector 
0 ER | 


The firm meets the acceptance criteria. We estimate its 95 % lower confidence bound 


to be 


aT 


LCB = 0.8222825 . 


Because LCB < LCBFB, we do not give a bonus to the firm. 
® Case # 3 


At the end of inspection, firm has following failure vector 
C20 Ol “OTF 


The firm meets the acceptance criteria. We estimate its 95 % lower confidence bound 


to be 


LCB = 0.8766983 . 


Because LCB > LCBFB, we give a bonus to the firm. 
2. Double Sampling Bonus System 
® Case # 1 


At the end of first inspection, firm has following failure vector. 


Cl Omori 


The firm meets the acceptance criteria. We estimate its 95 % lower confidence bound 


to be 
COBiy= 0822282502. 


Because LCB1 < LCBFB , we give a second chance to the firm. Here is the result of 


second inspection 


(H0m 3 lie o 


The aggregated failure vector will be 
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(ime 200 


It is clear that LCB2 < LCBI for this data, thus we do not give a bonus to the firm. 
sw Case # 2 


At the end of first inspection, firm has following failure vector 
fb 0-20. 1 ). 


The firm meets the acceptance criteria. We estimate its 95 % lower confidence bound 


to be 
LCB1 = 0.8222825 . 


Here, LCB < LCBFB thus we will give a second chance to the firm. Here is the result 


of second inspection 


(0: 00 tet); 


The aggregated failure vector will be 


and we estimate its 95 % lower confidence bound to be 
LCB1 = 0.8369564 . 


Because LCB < LCBFB, we do not give a bonus to the firm. 
= Case # 3 


At the end of first inspection, firm has following failure vector 


(Ch 00. 4), 
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The firm meets the acceptance criteria. We estimate its 95 % lower confidence bound 


to be 


LCB! = 0.822282) . 


Because LCB|1 < LCBFB, we will give a second chance to the firm. Here is the result 


of second inspection 


(0.020 


The aggregated failure vector will be 


G0 Ol 


and we estimate its 95 % lower confidence bound to be 
LCB1 = 0.8706521 . 


Because LCB2 > LCBFB, we will give a bonus to the firm. 


60 


VI. SIMULATION RESULTS OF BONUS SYSTEMS 


A. BACK GROUND 

We simulated the bonus system to get an idea of how well the bonus system pro- 
posed in the previous chapter works. This provides the user with a means of setting the 
cut off criteria for awarding a bonus. To generate random failure vectors, we need to 
know the probabilities of being successful in each test for the firm. First, we assume that 
the firm has an equal probability of being successful in each test, with the following 


Pero, — 1, — Pf, = P, 
e 0.9375 
e 0.9500 
¢ 0.9750 
e 0.9900 
e 0.9950. 


In a second set of simulations, we assume that the manufacturer test has a bigger 
probability of being successful than the other environment tests. The following proba- 


bilities of being successful in manufacturer test were used. 
e 0.9990 
e 0.9950 
¢ 0.9750 


e 0.9500. 
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In these cases, we used equal probabilities of being successful in other environment tests 
with values between 0.9375 and manufacturer test probability for that given case. 
Finally, we assume that the manufacturer test should have smaller probability of being 
successful. When we study these cases with this assumption, we used probabilities of 
being successful in other environment tests between manufacturer test probability and 
0.9990 for that given case. We assume that minimum value of the manufacturer test 
probability will be 0.9500 because of worst case in last two assumptions. 

We generated 2000 random failure vectors for each possible combination of proba- 
bilities by using program RANVEC in Appendix B. They were used for 1000 repli- 
cations of each bonus system, because the Double Sampling Bonus System ( DSBS ) 
potentially uses two failure vectors per replication. 

After getting the failure vectors, the next step is to decide lower confidence bound 
for giving bonus ( LCBFB ). We chose 0.800, 0.825, 0.850, 0.875, 0.900, 0.950, 0.999 
as LCBFB during our simulations. We simulated bonus systems by using program BO- 
NUS in Appendix G. Program BONUS counts how many times firm gets bonus during 
1000 replications. And then it calculates the bonus percentage dividing counted number 
by 1000. For each scenario this bonus percent is an estimate of the probability of getting 
a bonus. 

B. INITIAL COMPARISON OF SYSTEMS 
In this part, with equal probabilities of being successful in tests, we tried to see the 


difference between the Single and the Double Sampling Bonus Systems. For this reason, 


we simulated Single Sampling Bonus System. Results were tabulated and plotted in next 


pages. 


62 


Table 12. SINGLE SAMPLING BONUS SYSTEM (EQUAL PROBABILITIES) 


PROB.’S CHOSEN LOWER CONFIDENCE BOUNDS FOR BONUS 


[ome [os [oe [ons [ove [se [oom 
Toms [oom [ems [eos [em [one [ane [eon 
Toss [oom [eo ows [ous [we [eo [oor 
eo 
[sve [om [esi esr [aor [woe [oo [om 


It is obvious that there is no difference between some lower confidence bounds for bonus 


from the table above. The firm gets the same bonus percentage when we use 0.825, 
0.850, 0.875 as lower confidence bound for bonus. The same thing occurs when we use 
0.900 and 0.950 as LCBFB. For this reason, we are going to see four curves in Figure 


1. Figure 2 shows Double Sampling Bonus System with different LCBFB’s. 
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SINGLE SAMPLING BONUS SYSTEM 
DIFFERENT LCB'S FOR BONUS WITH EQUAL TEST PROBABILITIES 


ul 
5) 
e 
wi 
a 
vi 
2 
z 
o 
a 


Figure 1. Single Sampling Bonus System With Different LCBFB’s 


DOUBLE SAMPLING BONUS SYSTEM 
OIFFERENT LCB'S FOR BONUS WITH EQUAL TEST PROBABILITIES 


BONUS PERCENT 


Figure 2. Double Sampling Bonus System With Different LCBFB’s 


As you see in figure 2, we have six different curves for Double Sampling Bonus System. 
Because in this system, we can see from the table that LCBFB’s 0.900 and 0.950 have 


approximately the same bonus percentage. For this reason we did not plot for 0.950. 


Table 13. DOUBLE SAMPLING BONUS SYSTEM (EQUAL PROBABILITIES) 


PROB.’S CHOSEN LOWER CONFIDENCE BOUNDS FOR BONUS 


ome [ons [os [ns [a [oe om 
Tavs [vme [oom [ons [om [ene [oom [oom 
Tasso [oo [ome [ ons [we [oom [ com 
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sve [ore [ome [ee [oer [ow [wan [om 
Tere [oom [ome [orm [om [orm [os [nar 


Figure 2 shows DSBS to be more sensitive in the sense that a higher percentage of 


the firms were awarded a bonus. For this reason, we decided to implement DSBS in all 


simulations. 
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C. SIMULATION RESULTS WITH DIFFERENT LCB’S FOR BONUS 
In this section, results for each chosen lower confidence bound for bonus will be 
presented as follows. We used values in Table 13 to draw plots with equal probabilities 


in each test. Polynomial approximation was used in curve fitting. 
DOUBLE SAMPLING BONUS SYSTEM 
EQUAL TEST PROBABILITIES; LCBFB = 0.800 
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Figure 3. Double Sampling Bonus System With LCBFB = 0.800 


Table 14. DSBS (EQUAL PROBABILITIES) LCBFB = 0.800 


LOWER CONFIDENCE BOUND FOR BONUS IS 0.800 


[tow cosnpesce noun FoRoNs Sem 
ame [oe [eve [om | 
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In Figure 3, when the probabilities of being successful in each test increase, then the 
bonus percent increases. For example, when the test probabilities is equal to 0.9200, 
then the bonus percent is 0.026. If the firm increases its probabilities of being succesful 
in each test to 0.9900, then the bonus percent becomes 0.948. The result of the double 
sampling bonus system with different lower confidence bounds are tabulated and plotted 
in Appendix J. 

To see whether the probability of passing the manufacturer test effects the bonus 
percentage differently than the other probabilities, we assumed that the environment 
tests would have equal probabilities. But the firm will have a different probability of 
being successful in the manufacturer test. In the following table, the first row represents 
probability of being successful in the environment tests and the first column represents 
probability of being successful in the manufacturer test. The intersection of rows and 
columns gives us the bonus percentage of a firm with given probabilities of being suc- 
cessful in the tests. These procedures were done for each LCBFB value separately. The 
results are tabulated and plotted in Appendix K. Bonus percentages for LCBFB = 
0.8000 are plotted on the following page. It is clear that the bonus percentage of the 
firm will be high if the firm has big probabilities of being succesful in both manufacturer 


test and other joint environment tests. 
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Bonus percentages are tabulated and ploted with different probabilities. 


LCBFB = 0.800 


DOUBLE SAMPLING BONUS SYSTEM 
LCBFB = 0.800 
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Figure 4. Double Sampling Bonus System With LCBFB = 0.800 


Table 15. DSBS (DIFFERENT PROBABILITIES) LCBFB = 0.800 
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D. BONUS PERCENTAGE (BPRCT ) FORMULATION 

In this section we try to approximate probability of getting a bonus as a function 
of lower confidence bound for bonus ( LCBFB ), probability of passing manufacturer 
test ( P, ) , and probability of passing environment tests ( P, ). To do this, we use re- 
gression analysis with GRAFSTAT using the bonus percentages which are the simulated 
values of the probability of getting a bonus. Ifa reasonable relationship between the 
bonus percentage and LCBFB and the probabilities P, and P, can be found; it can be 
used to set LCBFB without resorting to simulation. 

We want to formulate bonus percentages as a function of the following: 

= Lower confidence bound for bonus ( LCBFB ) 

= Probability of passing manufacturer test P, 

s Probability of passing environment tests P,. 

After polynomial approximation, we can see that plotted graphs (Figures 3-10) look 
like logistic growth curves [Ref. 12: p. 383]. Then we can formulate bonus probabilities 


(BP) using the logistic growth function. 
l 
BP = ———— (7.4) 
where 
A = By + By Py + Bo Py + By Py + By Pi + Bs Py Py + Bo LCBFB (7.5) 


We can make a transformation described as below. 


l+e4= a. (7.6) 
fee 1 (7.7) 
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aoalos | boat I (7.8) 


Now we have a linear equation as a function of LCBFB, P, and P,. We may use 
simulation results and approximate parameters doing a linear regression. 

« Let Y be the (n x 1) column vector of observations on dependent variable. 

= Let X be the (nx p’ ) matrix consisting of a column of ones, which is labeled 1, 
followed by p column vectors of observations on independent variables. 

= Let B be the ( p’ x 1 ) vector of parameters to be estimated. 


= Let ¢ be the (nx 1 ) vector of random errors. Then 
Y= AX 0-22 (7.9) 
We can obtain observations on dependent variable as below, 


1 — BPRCT 
Y= log | =P DRC Tae | (7.10) 
where BPRCT is the bonus percentage obtained from the simulation. As an example 


we have the following information from simulations. 


« LCBFB = 0.875 


- P, = 0.9500 
- P? = 0.9025 
- P, = 0.9375 
- P? = 0.8789 


- P, P, = 0.8910 
- BPRCT = 0.029 
~Y = 3.5110. 


We can write following equation. 
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3.5110 = By + B, 0.95 + B, 0.9025 + B, 0.9375 + B, 0.8789 + B; 0.8910 + f.0.875 + €. 


We deliberately chose 30 random results from simulations. We formulate them as 
above. We can estimate unknown parameters on these equations doing linear re- 


gression. 30 equations can be written with matrix notation as follows: 


¥(1) 
¥(2) 
= | (7.11) 


¥(30) 
And X matrix is designed as follows: 


i Pay Pil) Pal) Ps(l) P, P,(1) 
1 Pi(2) P2(2) Po(2) P3(2) Py P2(2) 
X= | 1 PG) Pi) PGB) P33) Py P23) |. (7.12) 


1 P,(30) P7(30) P(30) P3(30) P, P,(30) 
And parameter vector will be as below. 


Bo 
B, 
By 
=| Bs. (7.13) 
Ba 
Bs 
Be 
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We estimate parameters with linear regression using GRAFSTAT packages. And 
then we formulate the bonus percentage using the estimates of parameters. The esti- 


mated bonus probabilities are given by 


Pee a Lae (7.14) 
1+ e%8 


The following estimates of the parameters were obtained from the regression. 


Table 16. COEFFICIENT VECTORS OF REGRESSION ANALYSIS 


Lo 2 


The standard error was 0.45471 after the linear regression. We can use these esti- 


mates of parameters in equation (7.5). And we can estimate the bonus probability of the 
firm when the probability of passing tests and lower confidence for bonus are known. 
Some examples are calculated with equation (7.5). Results of the calculations and 


comparison with simulation are summarized in following table. 


v2 


Table 17. BONUS PERCENTAGES WITH REGRESSION ANALYSIS 


ce 
0.800 0. | 0.9750 | 0. | 09750 | 0. | oes | 0. | 0564 
Se 
a 
a 
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VII. CONCLUSIONS AND RECOMMENDATIONS 


We estimated the reliability of pyrotechnic device from the sampling plan using the 
following models: 

« Maximum likelihood assuming independence of pyrotechnic device activation in 
different environments; 

« Log linear model incorparating some dependence between pyrotechnic device ac- 
tivation in different environments; 

« Worst case scenario which gives a lower bound for the estimated reliability for the 
general model where no assumptions are made about the form of dependence between 
pyrotechnic device activation in different environments. 

Using these models and based on sampling plan data, estimates of overall reliability 
along with 95 % lower confidence bound are obtained. We computed the lower confi- 
dence bound for each possible case by bootstrapping. 

Results from the first model are not consistent with the way pyrotechnic devices 
Operate. In particular, the estimated reliabilities are not ordered as we expected them 
to be. Intiutively, we expect samples with fewer failures to give smaller R values and 
corresponding LCB’s than samples with more failures. This is not the case for the first 
model. For example, the failure vector that has the maximum number of failures for 
each test ( a total of 5 failures ) is in the middle of the order with respect to R and the 
95 % lower confidence bound. The failure vector (0 1 0 1 ) with a total of 2 failures 
has a much lower R and LCB. The discrepency between the results of this model and 
what we expected to see are probably due to the fact that there is a dependence between 


the events which a device activates under different environments. 
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The results of the log linear model design with two-way dependence assumptions are 
similar to the first model. Two-way interaction terms Were used in this model. More- 
over, reliabilities of several cases were not calculated because of computational limita- 
tions. Due to the extensive censoring in the sample data, it does not appear to be 
possible to estimate the reliability based on models which incorparate dependence. 
Thus, we turn to finding lower bounds for the estimated reliability based on models with 
dependence. 

Finally, the worst case scenario model gives the most reasonable results for this 
problem. Both the lower bounds for estimates of reliabilities and 95 % lower confidence 
bounds for these lower bounds are ordered according to the total number of failures. 
Thus, the results from the worst case scenario were used to implement the bonus system 
for pyrotechnic devices. 

After getting the estimate of reliability and 95 % lower confidence bound for each 
case, we tried to design a bonus system to improve the quality of pyrotechnic devices. 
We used 95 % lower confidence bounds instead of the estimated overall reliabilities to 
decide whether to give bonuses. Two sampling plans for giving the bonus to the man- 
ufacturers were considered: 

# Single Sampling Bonus System 

®# Double Sampling Bonus System. 

We simulated two sampling bonus systems to see the difference between them. 
We concluded that the double sampling system is more sensitive than single sampling 
bonus system. 

To formulate an approximate bonus percentage as a function of lower confidence 
bound for bonus ( LCBFB ), the probability of passing manufacturer test P,, and the 
probability of passing the environment tests P,, we used regression analysis with 


GRAFSTAT. 
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Because the sampling plan results in so much censoring, the only reasonable esti- 
mates of overall reliability that we obtained were actually lower bounds. This makes the 
bonus system conservative in the sense that a bonus might not be awarded what it is 
deserved. Therefore, if a bonus system is to be implemented, a more comprehensive 
sampling plan needs to be devised which allows estimation of R. A simple solution to 
this problem can be to apply all environmental tests to the same sample which would 


give a measure of dependence between environmental tests. 
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APPENDIX A. PROGRAM MLEA 


RRERRERRRRRREEREERRREERRREEREREEEREREREREREREERRERRRERERERERERERRREEREREE 


PROGRAM MLEA 
HRARKRARERKRAERERERERE RRR EER ERERERER ER RRERERE RRERERERRRERERERRRERE 

THIS IS A FORTRAN PROGRAM TO CALCULATE THE RELIABILITY OF AN 
ITEM AFTER EXPOSURE TO SEVERAL ENVIRONMENTS WITH INDEPENDENCE 
ASSUMPTION WHICH IS DESCRIBED IN CHAPTER I . THE PROGRAM READS 
5000 SUCCESS VECTORS, WHICH ARE RANDOMLY GENERATED BY PROGRAM 
RANVEC, FROM AN INPUT FILE CALLED SUCVECT ONE BY ONE. AFTER 
CALCULATION IT THEN PROGRAM WRITES ESTIMATED RELIABILITIES TO 
AN OUTPUT FILE CALLED RESULT. 


FOIE AL IRIE AR EL EER IL IL IL AL IL AE IL AE IL IE GEIR ERIE IN IL IR IL EE IN GD INGER GEIR GL IL ITED IT ITIL IVER IL ER ILIV GRID IL AIL ITIL IL ITIL ED IL GL ID ID GD OD 


VARIABLES 

witkkickiickiiK CIARA RAIA IRI RARE RRR ICI 

SOM : NUMBER OF SUCCESSFUL ITEMS IN MANUFACTURER TEST. 

SOTH : NUMBER OF SUCCESFUL ITEMS IN TEMPERATURE AND HUMIDITY 
TEST. 

SOV : NUMBER OF SUCCESSFUL ITEMS IN VIBRATION TEST. 

SOA : NUMBER OF SUCCESSFUL ITEMS IN ALTITUDE TEST. 

R1H : ESTIMATED PROBABILITY OF PASSING FROM MANUFACTURER TEST. 

R2H : ESTIMATED PROBABILITY OF PASSING FROM TEMPERATURE AND 
HUMEDITY*TEST . 

R3H : ESTIMATED PROBABILITY OF PASSING FROM VIBRATION TEST 

R4H : ESTIMATED PROBABILITY OF PASSING FROM ALTITUDE TEST 


RHMLE : ESTIMATED RELIABILITY OF ITEM AFTER EXPOSURE TO 
SEVERAL ENVIRONMENT TESTS. 

N : SAMPLE SIZES FOR EACH TEST 

X : NUMBER OF SUCCESSFUL ITEMS IN EACH TEST. 

FLAG : INDICATOR VARIABLE FOR DETERMINING EASY AND HARD CASE. 

R1IMAX : R1H VALUE WHICH MAXIMIZES LIKELIHOOD. 
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R2MAX : R2H VALUE WHICH MAXIMIZES LIKELIHOOD. 
R3MAX : R3H VALUE WHICH MAXIMIZES LIKELIHOOD. 
R4MAX : R4H VALUE WHICH MAXIMIZES LIKELIHOOD. 
LMAX : HIGHEST MAXIMUM LIKELIHOOD VALUE. 


L : LIKELIHOOD VALUES AT THE END OF EACH HARD CASE. 

V : TOTAL SAMPLE SIZES IN EACH HARD CASE. 

Q1H : REORGANIZED PROBABILITY OF MANUFACTURER TEST. 

Q2H : REORGANIZED PROBABILITY OF TEMP. AND HUMIDITY TEST. 
Q3H : REORGANIZED PROBABILITY OF VIBRATION TEST. 

Q4H : REORGANIZED PROBABILITY OF ALTITUDE TEST. 


keke eek bebe ERA REAR RARE ERR REAR ER ER ER EEA AER ER RERER ER 


TYPE DECLARATION 
REAL SOM(5000) ,SOTH(5000) ,SOV(5000) ,SOA(5000) ,R1H,R2H,R3H,R4H, 
+RHMLE(5000) ,N(4) ,X(4) , RIMAX, R2MAX , R3MAX, R4MAX , LMAX,V(7) ,Q1H,Q2H, 
+Q3H ,Q4H ,NTOT,L(7) 

INTEGER I,J 


PRR EERE ICT ICIS TOTTI TN ISIS TS ICTS TS OTS TS TS TTS TE TST TC TE TS TOBE ETS TS STE TS ICTS TS TS TOBE 


READING SUCCESS VECTORS FROM SUSVECT FILE 
DO 60 I=1,5000 
READ(7,*) SOM(I),SOTH (I),SOV(I),SOA(I) 


SITTIN TS FSGS OTS TCT TS 


eae 


FILES FOR READING AND WRITING 
CALL EXCMS (*FILEDEF 7 DISK SUCVECT DATA Al") 
CALL EXCMS ('FILEDEF 16 DISK RESULT DATA Al") 


INITIALIZATION OF SAMPLE SIZES 
N(1)= 20.0 
N(2)= 20.0 
N(3)= 32.0 
N(4)= 20.0 
INITIALIZATION OF FLAG VARIABLES 
DO 10 J=2,4 
FLAG(J)= .FALSE. 


78 


10 CONTINUE 
NUMBER OF SUCCESSES IN EACH TEST 
X(1)= SOM(I) 
X¥(2)= SOTHC) 
X(3)= SOV(I) 
X(4)= SOACI) 
Setedetetededetedotedete dete RREEREREEEREEREREEERRRRRERERERERRRERRERRERREAE 
CHECK OPERATION FOR EASY CASE 
IF ((X(1)/20.00).GE.(X(2)/20.00)) THEN 
FLAG(2)= .TRUE. 
END IF 
IF ((X(1)/20.00) .GE.(X(3)/32.00)) THEN 
FLAG(3)= .TRUE. 
END IF 
IF ((X(1)/20.00).GE.(X(4)/20.00)) THEN 
FLAG(4)= .TRUE. 
END IF 
IF ( FLAG(2).AND.FLAG(3).AND.FLAG(4) ) THEN 
tedetedoddcictedede RR REAR RREREERRREERERRRERREREERRECRERREERRREEERERRERE 
CALCULATIONS IN EASY CASE 
R1H= X(1)/20.0 
R2H= X(2)/20.0 
R3H= X(3)/32.0 
R4H= X(4)/20.0 
Q1H= R1H 
Q2H= R2H/Q1H 
Q3H= R3H/Q1H 
Q4H= R4H/Q1H 
RHMLE(I)= Q1H*Q2H*Q3H*Q4H 


GO TO 50 


CALCULATIONS IN HARD CASES 
INITIALIZATION 
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30 


40 


DO 30 J=1,7 
V(J)= 0.0 
L(J)= 0.0 
CONTINUE 
LMAX= 0.0 
NTOT= 0.0 
Heekieteddek tek RARER ARERERRERRRE ERR ER ERR ERR ARERR REERRERRERREREREREERE 
CASE 1 R1H.LT.1.0 AND R1H = R2H = R3H = R4H 
DO 40 J=1,4 
V(1)= V(1)+X(J) 
NTOT= NTOT+N(J) 
CONTINUE 
R1H= V(1)/NTOT 
L(1)= ((R1H)**V(1))*((NTOT-V(1)) /NTOT)**(NTOT-V(1)) 
IF (L(1).GT.LMAX) THEN 
LMAX= L(1) 
R1MAX= R1H 
R2MAX= R1H 
R3MAX= R1H 
R4MAX= R1H 
END IF 


@eevVevevavavayayae 


CASE 2 R1H.LT.1.0, R1H = R2H = R3H , R4H IS BETWEEN 0.0 AND R1H 
IF ( X(4)/20.0.NE.1.0 ) THEN 
V(2)= X(1)+X(2)+X(3) 
NTOT= N(1)+N(2)+N(3) 
R1H= V(2)/NTOT 
IF((X(4)/20.00).LE.R1H) THEN 
L(2)= ((R1H)**V(2))*( (NTOT-V(2)) /NTOT)**(NTOT-V(2) )*((X(4) 


+/20.0)**X(4))*(((20.0-X(4) )/20.0)**(20.0-X(4))) 


IF (L(2).GT.LMAX) THEN 
LMAX= L(2) 
R1MAX= R1H 
R2MAX= R1H 
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R3MAX= R1H 
R4MAX= X(4)/20.0 
END IF 
END IF 
END IF 


Here dee RARER RR ARRER ARR ERE RRR RRERR ER RRERRRERERERE RE ERERR ER REERE 
CASE 3 R1H.LT.1.0, R1H = R2H = R4H , R3H IS BETWEEN 0.0 AND R1H 


IF ( X(3)/32.0.NE.1.0 ) THEN 
V(3)= X(1)+X(2)+X(4) 
NTOT= N(1)+N(2)+N(4) 
R1H= V(3)/NTOT 
IF ((X(3)/32.00).LE.R1H) THEN 
L(3)= ((R1H)**V(3))*((NTOT-V(3) )/NTOT)**(NTOT-V(3) )*((X(3) 
+/32.0)**X(3))*(((32.0-X(3))/32.0)**(32.0-X(3))) 
IF (L(3).GT.LMAX) THEN 
LMAX= L(3) 
R1MAX= R1H 
R2MAX= R1H 
R3MAX= X(3)/32.0 
R4MAX= R1H 
END IF 
END IF 


END IF 
RRERRERERRERERREREREREEERERERERERERERERERERERERERERERRERERERARRREER 


CASE 4 R1H.LT.1.0, R1H = R3H = R4H , R2H IS BETWEEN 0.0 AND R1H 
IF ( X(2)/20.0.NE.1.0 ) THEN 
V(4)= X(1)+X(3)+X(4) 
NTOT= N(1)+N(3)+N(4) 
R1H= V(4)/NTOT 
IF ((X(2)/20.00).LE.R1H) THEN 
L(4)= (C(R1H)**V(4) )*(C(NTOT-V(4) ) /NTOT)**(NTOT-V(4) )*((X(2) 
+/20.0)**X(2))*(((20.0-X(2)) /20.0)**(20.0-X(2))) 
IF (L(4).GT.LMAX) THEN 
LMAX= L(4) 
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R1IMAX= R1H 
R2MAX= X(2)/20.0 
R3MAX= R1H 
R4MAX= R1H 
END IF 
END IF 
END IF 
* CASE 5 R1H.LT.1.0, R1H = R2H R3H , R4H ARE BETWEEN 0.0 AND R1H 
IF ((X(3)/32.0.NE.1.0).AND.(X(4)/20.0.NE.1.0)) THEN 
V(5)= X(1)+X(2) 
NTOT= N(1)+N(2) 
R1H= V(5)/NTOT 
IF((X(3)/32.00) .LE.R1H.AND.(X(4)/20.00).LE.R1H) THEN 
L(5)= ((R1H)**V(5))*((NTOT-V(5) )/NTOT)**(NTOT-V(5))*((X(3) 
+/32.0)**X(3))*((32.0-X(3))/32.0)**(32.0-X(3))*((X(4)/20.0)** 
+x(4))*((20.0-X(4) )/20.0)**(20.0-X(4)) 
IF (L(5).GT.LMAX) THEN 
LMAX= L(5) 
R1MAX= R1H 
R2MAX= R1H 
R3MAX= X(3)/32.0 
R4MAX= X(4)/20.0 


END IF 
END IF 
END IF 
* FARRAR AEE ERA IARI IIIA REI RRR REE RRERERE EE 
* CASE 6 R1H.LT.1.0, R1H = R3H R2H , R4H ARE BETWEEN 0.0 AND R1H 


IF ((X(2)/20.0.NE.1.0).AND.(X(4)/20.0.NE.1.0)) THEN 
V(6)= X(1)+X(3) 
NTOT= N(1)+N(3) 
R1H= V(6)/NTOT 
IF((X(2)/20.00).LE.R1H.AND.(X(4)/20.00).LE.R1H) THEN 
L(6)= ((R1H)**V(6))*((NTOT-V(6) ) /NTOT)**(NTOT-V(6) )*((X(2) 
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++ 


+/20.0)**X(2))*((20.0-X(2))/20.0)**(20.0-X(2))*((X(4)/20.0)** 
+x(4))*((20.0-X(4))/20.0)**(20.0-X(4)) 


IF (L(6).GT.LMAX) THEN 
LMAX= L(6) 
R1MAX= R1H 
R2MAX= X(2)/20.0 
R3MAX= R1H 
R4MAX= X(4)/20.0 
END IF 
END IF 
END IF 
KEKKANRERARREKR EERE RR ERRERER REE EER ERE RRR EER EREERERREERERERER ERE 
CASE 7 R1H.LT.1.0, R1H = R4H R2H , R3H ARE BETWEEN 0.0 AND R1H 
TE GOcC2) / 20. 02NE. 1.0) AND? (X03) /32.0.NE.1.0)) THEN 
V(7)= X(1)+X(4) 
NTOT= N(1)+N(4) 
R1H= V(7)/NTOT 
IF((X(2)/20.00) .LE.R1H.AND.(X(3)/32.00).LE.R1H) THEN 
L(7)= ((R1H)**V(7))*CCNTOT-V(7) ) /NTOT)**(NTOT-V(7))*((X(3) 
ep o230) 523) )*((3200-XC3) )/32.0)**(32.0-X(3))*CCX(2) 720.0) =*k(@2)) 
tt e2 0). O-X( 27 20.0) "(20 .0-X02)) 


IF (L(7).GT.LMAX) THEN 
LMAX= L(7) 
R1IMAX= R1H 
R2MAX= X(2)/20.0 
R3MAX= X(3)/32.0 
R4MAX= R1H 

END IF 

END TF 

END IF 


whealsaloale alrals ale al, alaagkei 
RRB BATA Ie SNES ENTS TS DCTS ISIS TS TS TS TS TS TS TS TS OTT TE TTT OER ERERERERRERERERE 
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50 
60 


CALCULATION OF RHMLE VALUE WITH RESPECT TO CASE WHICH HAS LARGEST 


MAXIMUM LIKELIHOOD IN HARD CASE 


R1H= 
R2H= 
R3H= 
R4H= 
Q1H= 
Q2H= 
Q3H= 
Q4H= 


R1MAX 
R2MAX 
R3MAX 
R4MAX 
R1H 
R2H/Q1H 
R3H/Q1H 
R4H/Q1H 


RHMLE(I)= Q1H*Q2H*Q3H*Q4H 


REECE EREEREREEREREEREEREREREREERREREEREERERRRERERERRERE 


WRITING AFTER EACH CALCULATION 
WRITE (16,*) RHMLE(1) 
CONTINUE 


STOP 
END 
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APPENDIX B. PROGRAM RANVEC 


hake abe aks leaks abs ale ale akeakentaatealeals Lealents abouts aleaks ake als oleate abeats alealealents afeakes abe tents aloalsals alealsale ale ale aks leaks ate mts ks abs ale abs nto aleals bealselseliatonts Le 
FUER OV ON ONAN EN ENOL EVER IRON ER GL IN IN FRET GL IN GC ER TR ECTC ERECT IN EN FUER OD IU ID TIRED TLIC EVIL IN IN IL GN IN ER TV EV EC PCN IN ED ON ERIN ONIN IN GSU OD 


PROGRAM RANVEC 

THIS IS THE PROGRAM TO GENARETE RANDOM NUMBERS FROM BINOMIAL 
DISTRIBUTION.THE PROGRAM READS PROBABILITIES OF BEING SUCCESSFUL 
IN FOUR TESTS INTERACTIVELY. IT GENERATES UNIFORMLY DISTRIBUTED 
RANDOM NUMBERS WITH THESE PROBABILITIES ACORDING TO SAMPLE SIZE 
OF EACH TEST. FOR EACH TEST PROGRAM COUNTS UNIFORMLY DISTRIBUTED 
RANDOM NUMBERS, WHICH HAVE GREATER THAN OR EQUAL PROBABILITY 
WITH RESPECT TO THE GIVEN PROBABILITY FOR THAT TEST.TOTAL COUNTS 
GIVE US SUCCESSFUL ITEM NUMBERS FOR EACH TEST. PROGRAM UPDATES 
SEEDS AND CALLS SUBROUTINE RANNUM IN EACH ITERATION. THE PROGRAM 
GENERATES 5000 SUCCESS VECTOR AND WRITES THEM TO AN OUTPUT FILE 
CALLED SUSVECT. 


ple EP et Peat Pe at Seat Ba et et nl wlevloulentealoate tects tectoctontoaRecBsstectoats ckeetaclontantontertan tents 
WNBA TNS STIR TST TORR TTR TER RRR ARERRRBRERBRRRRRARRRRRRERRREEREEREERRER 


VARIABLES .- 

PIKERICR RIE ER IRIE ERE REE ERE RII ER ER EER ER ERE 

PSIM : PROBABILITY OF SUCCESS IN MANUFACTURER TEST. 

PSITH : PROBABILITY OF SUCCESS IN TEMP. AND HUMIDITY TEST. 

Poly : PROBABILITY OF SUCCESS IN VIBRATION TEST 

PSIA : PROBABILITY OF SUCCESS IN ALTITUDE TEST 

NUM1 : COUNTER FOR MANUFACTURER TEST 

NUM2 : COUNTER FOR TEMP. AND HUMIDITY TEST 

NUM3 : COUNTER FOR VIBRATION TEST 

NUM4 : COUNTER FOR ALTIUTDE TEST 

X : BINOMIALLY DISTRIBUTED RANDOM NUMBER FOR MANUFACTURER 
TEST 

V : BINOMIALLY DISTRIBUTED RANDOM NUMBER FOR TEMP.AND 


HUMUDITY TEST 
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Y : BINOMIALLY DISTRIBUTED RANDOM NUMBER FOR VIBRATION TEST 
: BINOMIALLY DISTRIBUTED RANDOM NUMBER FOR ALTITUDE TEST 


A : UNIFORMLY DISTRIBUTED RANDOM NUMBER FOR MANUFACTURER 
TEST 

B : UNIFORMLY DISTRIBUTED RANDOM NUMBER FOR TEMP.AND 
HUMIDITY TEST 

C : UNIFORMLY DISTRIBUTED RANDOM NUMBER FOR VIBRATION TEST 

D : UNIFORMLY DISTRIBUTED RANDOM NUMBER FOR ALTITUDE TEST 


ISEED : SEED NUMBER FOR MANUFACTURER TEST 

KSEED : SEED NUMBER FOR TEMP. AND HUMIDITY TEST 

LSEED ;: SEED NUMBER FOR VIBRATION TEST 

MSEED : SEED NUMBER FOR ALTITUDE TEST 

Sevedede tee RRR RRR ERRRRERERRRRRRRE RRR ERRRRRERERERIRERE 
TYPE DECLARATION 

REAL PSIM,PSITH,PSIV,PSIA,NUM1,NUM2,NUM3 ,NUM4 ,X(5000) ,V(5000) 
+,Y(5000) ,Z2(5000) ,A,B,C,D 

INTEGER ISEED,KSEED,LSEED ,MSEED 

PRR IPR RRR ICTR REAPER RRR ERAR TERRIER 
INITIALIZATION 

ISEED = 45267 

KSEED = 113234 

LSEED = 435 

MSEED = 1 
PARKER RIKER EREREEEEEREREERERERREEREER 
READING TEST PROBABILITIES OF BEING SUCCESSFUL IN EACH TEST 
WRITE(*,*)'PLEASE WRITE PROBABILITY OF BEING SUCCESSFUL IN 
+ MANUFACTURER TEST' 

READ (*,*) PSIM 

WRITE(*,*)'PLEASE WRITE PROBABILITY OF BEING SUCCESSFUL IN 
+ TEMPERATURE AND HUMIDITY TEST' 

READ (*,*) PSITH 

WRITE(*,*)'PLEASE WRITE PROBABILITY OF BEING SUCCESSFUL IN 
+ VIBRATION TEST' 

READ (*,*) PSIV 
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WRITE(*,*)'PLEASE WRITE PROBABILITY OF BEING SUCCESSFUL IN 

+ ALTITIDE TEST 

READ (*,*) PSIA 

RRR RRERERERRERERERR ERA RRETERRRRRRERRRERRRERERERRRRERRRRRRRER: 

OPENING AN OUTPUT FILE TO WRITE THE RESULTS 

CALL EXCMS ('FILEDEF 16 DISK SUSVECT DATA Al1') 

Ree KERR RRR RRR RERRRREERERRRRERERERREREREERR RR RR 

GENERATION 

DO 50 J=1,5000,1 

INITIALIZATION IN EACH ITERATION 
NUM1= 0.0 
NUM2= 0.0 
NUM3= 0.0 
NUM4= 0.0 

PRRREKRRERERRERERRRAKERERERRERERRRERRRRRRRRERRERRERRRRRRERERRRRRERERRE 

SEEDS UPDATATION IN EACH ITERATION 
ISEED=ISEED+17 


KSEED=KSEED+1356 
LSEED=LSEED+1 
MSEED=MSEED+789 


GENERATION OF SUCCESSFUL NUMBER OF ITEMS IN MANUFACTURER TEST 
DO 10 I=1,20,1 
CALL RANNUM (1,ISEED,0.0,1.0,0.0,A) 
IF ( A.LT.PSIM ) THEN 
NUM1=NUM1+1 


END IF 
X(J) = NUM1 
CONTINUE 


RIKKI EEREEE REE ERR ERR ERRERE 
GENERATION OF SUCCESSFUL NUMBER OF ITEMS IN TEMP. AND HUM. TEST 
DO 20 K=1,20,1 
CALL RANNUM (1,KSEED,0.0,1.0,0.0,B) 
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ote 
rr 


30 


ae 
ry 


40 


50 


IF ( B.LT.PSITH ) THEN 
NUM2=NUM2+1 
END IF 
V(J) = NUM2 
CONTINUE 
GENERATION OF SUCCESSFUL NUMBER OF ITEMS IN VIBRATION TEST 
DO 30 Leleae 
CALL RANNUM (1,LSEED,0.0,1.0,0.0,C) 
IF ( C.LT.PSIV ) THEN 
NUM3=NUM3+1 
END IF 
Y(J) = NUM3 
CONTINUE 
PPR RERERERRERERRREERERRERREREERREREEERRRERRECESE 
GENERATION OF SUCCESSFUL NUMBER OF ITEMS IN ALTITUDE TEST 
DO 40 M=1,20,1 
CALL RANNUM (1,MSEED,0.0,1.0,0.0,D) 
IF ( D.LT.PSIA ) THEN 
NUM4=NUM4+1 
END IF 
Z(J) = NUM4 
CONTINUE 
DEPOT AED AIARRRERIEREREAEEEKERAEEEEREEREREEREREER 


WRITING THE RESULTS TO AN OUTPUT FILE AS 4 TUPLE 
WRITE (16,1) X(J),V(J),Y(J),2(J) 
FORMAT (1X,F12.7,4X,F12.7,6%,F12.7,4%)F12-7,.4%) 
CONTINUE 
STOP 
END 


SUBROUTINE RANNUM(DISTN, SEED, RPARM1, RPARM2, IPARM, X) 


WUT TS TTS TTT ICICI ICTS TS TODS TC TS ICTS TS TE ISTE TSS TST TS 


THIS SUBROUTINE IS A PART OF SIMUTIL FORTRAN WHICH IS WRITTEN 
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BY DR. M. P. BAILEY. THIS SUBROUTINE PROVIDES AN INTERFACE WITH 

THE LLRANDOMII ROUTINES PROVIDED IN THE NONIMSL LIBRARY. THE 

PARAMETER REQUIRMENTS AND CALLING PROCEDURES ARE AS FOLLOWS: 

DISTN = DISTRIBUTION TYPE YOU WANT TO SELECT AN INTEGER BETWEEN 1 
AND 7. 

SEED = THE RANDOM NUMBER SEED YOU WISH TO USE. 

RPARM1, RPARM2, AND IPARM ARE REAL AND INTEGER PARAMETERS PASSED 

TO THE ROUTINE WITH MEANINGS WHICH VARY WITH THE TYPE OF DISTRI_ 

BUTION YOU DESIRE. 

X = THE RETURNED RANDOM NUMBER, IT IS ALWAYS REAL. 

DISTRIBUTION NUMBERS AND THE ASSOCIATED PARM DEFINITIONS 

1--UNIFORM ON THE INTERVAL RPARM1 TO RPARM2. 

2--NORMAL WITH MEAN RPARM1 AND VARIANCE RPARM2. 

3--EXPONENTIAL WITH RATE RPARM1. 

4--COUCHY WITH A = RPARM1 AND B = RPARM2. 

5--GAMMA WITH SHAPE RPARM2 AND RATE RPARM1. 

6--POISSON WITH RATE RPARM1. 


7--GEOMETRIC WITH P = RPARM1. 
HAHEI ERARIR AEA ERRERAEREEREREERERRERRERRERE 


TYPE DECLARATION 
REAL RPARM1,RPARM2,X,TEMP,VARIAT(1) 


INTEGER DISTN, SEED, IPARM, N 
PIAA IIIA AIR AIA AIA RIERA ER ERE ER ERIE ERR EERIE 


IF (DISTN.LE.0.OR.DISTN.GT.8) THEN 
WRITE(10, *) ‘ILLEGAL CALL TO RANNUM, BAD DISTN' 
STOP 

ENDIF 

GOTONGIO, 20, 30, 40, 50, 60, 70), DISTN 


hiked? 
GENERATE A UNIFORM BETWEEN RPARM1 AND RPARM2 
CONTINUE 
IF (RPARM1 - RPARM2.EQ.0) THEN 
WRITE(10, *) ‘ILLEGAL EQUAL RPARMS IN RANNUM' 
STOP 
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20 


30 


ENDIF 
IF (RPARM1.GT.RPARM2) THEN 
TEMP = RPARM1 
RPARM1 = RPARM2 
RPARM2 = TEMP 
ENDIF 
CALL LRND(SEED, VARIAT, 1, 1, 0) 
VARIAT(1) = RPARM1 + (RPARM2 - RPARM1) * VARIAT(1) 


GOTO 80 


PEAR ARRREREREERERERREREREEREREERERRERERERRERRERRREREERREEREREREREEREERE 


GENERATE A NORMAL WITH MEAN RPARM1 AND STDDEV RPARM2 
CALL LNORM(SEED, VARIAT, 1, 1, 0) 

VARIAT(1) = (VARIAT(1) * RPARM2) + RPARM1 

GOTO 80 


on CXS MPO eT eT eee) vee ces 7 RR CT ee ee ey devs eae acids * cavieds PER AE TN ee Fe ae ed ah ae Ce ae ae CRTC COTO 


GENERATE AN EXPONENTIAL WITH RATE (1/MEAN) RPARM1 

CONTINUE 

IF (RPARM1.EQ.0) THEN 
WRITE(10, *) ‘ILLEGAL ZERO RATE IN RANNUM' 
STOP 

ENDIF 

CALL LEXPN(SEED, VARIAT, 1, 1, 0) 

VARIAT(1) = VARIAT(1) / RPARM1 


GOTO 80 


Rickiictes BIKER ERERREECEERERERREERERERERERERE 


GENERATE A COUCHY WITH A = RPARM1 AND B = RPARM2 

CONTINUE 

IF (RPARM2.LE.0) THEN 
WRITE(10, *) ‘ILLEGAL COUCHY SPREAD IN RANNUM, B = ',RPARM2 
STOP 

ENDIF 

CALL LCCHY(SEED, VARIAT, 1, 1, 0) 

VARIAT(1) = (VARIAT(1) * RPARM2) + RPARM1 
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80 


GOTO 80 
GENERATE GAMMA WITH SHAPE RPARM2 AND RATE RPRAM1 
CONTINUE 
IF (RPARM1.LE.0) THEN 
WRITE(10, *) ‘ILLEGAL NONPOSITIVE GAMMA RATE IN RANNUM' 
STOP 
ENDIF 
IF (RPARM2.LE.0) THEN 
WRITE(10, *) ‘ILLEGAL SHAPE PARAMETER IN RANNUM' 
STOP 
ENDIF 
CALL LGAMA(SEED, VARIAT, 1, 1, 0, RPARM2) 
VARIAT(1) = VARIAT(1) * (1.0 / RPARM1) 
GOTO 80 
FEST Ae ae OO OE OORT COE OE Tor oe eae SS OR ea oe ok oe eo ee ee a ce 
GENERATE POISSON WITH RATE RPRAM1 
CONTINUE 
IF (RPARM1.LE.0) THEN 
WRITE(10, *) ‘ILLEGAL POISSON RATE IN RANNUM' 
STOP 
ENDIF 
CALL LPOIS(SEED, VARIAT, 1, 1, 0, RPARM1) 
GOTO 80 
RREREREEREREREREREREERERAEERRERREERRERRRRERERERERERREREERERRERREEEE 
GENERATE GEOMETRIC WITH P = RPRAM1 
CONTINUE 
IF (RPARM1.LE.0) THEN 
WRITE(10, *) ‘ILLEGAL GEOM PROB IN RANNUM' 
STOP 
ENDIF 
CALL LGEOM(SEED, VARIAT, 1, 1, 0, RPARM1) 
GOTO 80 
CONTINUE 
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APPENDIX C. PROGRAM SORT 


ele bh Lebel ebeols leat ak label abealol Lialals heal als abeelials aleelalalials abeals lualealsals alsalsalsals aleals alsalratsats alsalsals aloals aloalsals ats ale leaks telat als 
FV GV EN TN GN GN EN GUN ON GE TL SLED GD ED GLEN TL EN GL ED ED EU TEIN YEN EN EN EN GY IVER ED ELEN EVID EVID IL IVIL ENT GLEN GN ID IU ENE EN TV IN GRIN ES ILIV GD ED OD 


PROGRAM SORT 

THIS IS THE SORTING PROGRAM. PROGRAM USES BUBBLE SORT ALGORITHM. 
PROGRAM READS ESTIMATED RELIABILITIES FROM AN INPUT FILE CALLED 
RESULT. IT SORTS FROM SMALLEST TO LARGEST, AND WRITES IN TO AN 
OUTPUT FILE CALLED FRESULT WITH 95 % LOWER CONFIDENCE BOUND. 


HERRERA ALAA AAAS A NANA AAE ERA R AR RER ERR R EERE RRR RR RRR ERE RR EERE RRR E RE 


VARIABLES 
Sesedededete ese RREREEREEEREREEEREEEREEREREREEEREREREREEREERE 
A : ESTIMATED RELIABILITY 
FLAG : INDICATOR VARIABLE TAKES VALUE ' OF ' AND ' OFF ' 
TYPE DECLARATION 
CHARACTER FLAG*3 
REAL A(5000) 
INTEGER I,N,J 
PRRREEEREREEERERERERREREREERREEEREREREEREREEREREREREREERERERREERERREEREERR 
OPENING AN INPUT AND AN OUTPUT FILE 
CALL EXCMS ('FILEDEF 9 DISK RESULT DATA A1') 
CALL EXCMS ('FILEDEF 15 DISK FRESULT DATA A1') 
EREREERRRERERERRERERERERRERRRRERERERRRRERRRREREREEERERREERRERERRRERREERRERE 
READING ESTIMATED RELIABILITIES 
DO 10 I=1,5000 
READ(9,*) A(I) 
CONTINUE 
SORTING OPERATION 
N=I-1 
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DOE30 1-N.2,-1 
FLAG='OFF' 
HO 20 J=15 i 
IF (A(J).GT.A(J+1)) THEN 


TEMP=A(J) 
A(J)=A(J+1) 
A(J+1)=TEMP 
FLAG='ON' 
END IF 
CONTINUE 
IF (FLAG.EQ.'OFF') THEN 
GO TO 40 
END IF 
CONTINUE 
CONTINUE 


= Fu eu ee ae eu ay 


WRITING THE RESULTS IN ASCENDING ORDER TO OUTPUT FILE 
DO 50 I=1,5000 
WRITE (15,*) ACI) 
CONTINUE 
WRITEM( 15:51) BEZ50) 
FORMAT (///,15X,'95 % LOWER CONFIDENCE BOUND IS',1X,F12.7) 
STOP 
END 
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APPENDIX D. PROGRAM INITIAL 


has hs hes ah ahs hes ahs abe ah ahs whe whe whe ahs ah ahs ahs whe ahs ae ahs ake aks whe ake ahs ake ake whe ake ales abe ake wks ahs ake ake ake aks ahoatc ake ale ake ahs ole ale whe ahs ahs ahs ahs ake aks ake aks whe ahs ake whe ake ake hs ake 
JNU AN GN IRAN ENGLER GR ILI IV IC GC ONGC IN ICING ES IRIN ISIN GL OCIA ID IL ERED GOEL OV EC EROS GC ED CIT IN GCI OSG IN ISIC IC IN PEIN EL TREC TS TN GS ES ES 


PROGRAM INITIAL 

THIS PROGRAM, CALCULATES INITIAL GUESSES FOR PARAMETERS IN 
LOGLINEAR MODEL BY MEANS OF PROGRAM PARAM WHICH IS IN APPENDIX E. 
PROGRAM SUPPLIES PARTIAL SUMS OF EXPECTATION, TO SOLVE EQUATIONS 
IN PROGRAM PARAM. IT READS INTERACTIVELY NUMBER OF FAILURES IN 
TESTS. PROGRAM WRITES RESULTS TO AN OUTPUT FILE CALLED EXPECT. 


whe he ake ahs ake hake ah ah aks ake whe ahs ahs ake ake ake ake ahewbo aks ake ake ake ake ahs ake whe ake ake ahs whe ake ake ahs ake abs ake ake ake ake howls ale mlc ale ahs ake uke ake Leake he ahs whe alente ahs ake ake wks ale ake aks 
FUG AN TE IN EN ON GR IN TS INGE PUGS IU INGUIN ON GD EU IN ILI ID IN IVIL GD IN IN TYING IN GY ICID IN IS IV ININ IL GU GV IN GN IV ISIN ILIV OLIN IDIVOV AN ET IL IN IN OD 


VARIABLES 

REAR RAERIRRIIRERAERAEARERER ERIE IRR ERIE RERRRER EERE 
FOM : NUMBER OF FAILURES IN MANUFACTURER TEST 

FOTH : NUMBER OF FAILURES IN TEMPERATURE AND HUMIDITY TEST 


FOV : NUMBER OF FAILURES IN VIBRATION TEST 

FOA : NUMBER OF FAILURES IN ALTITUDE TEST 

al : SUCCESS RATIO FOR MANUFACTURER TEST 

P2 : SUCCESS RATIO FOR TEMPERATURE AND HUMIDITY TEST 

RS : SUCCESS RATIO FOR VIBRATION TEST 

P4 : SUCCESS RATIO FOR ALTITUDE TEST 

Ql > FAILURE RATIO FOR MANUFACTURER TEST 

OZ : FAILURE RATIO FOR TEMPERATURE AND HUMIDITY TEST 

Q3 : FAILURE RATIO FOR VIBRATION TEST 

Q4 : FAILURE RATIO FOR ALTITUDE TEST 

X1 : PARTIAL SUM OF EXPECTED NUMBERS IN CELLS FOR EQUATION 1 
X2 : PARTIAL SUM OF EXPECTED NUMBERS IN CELLS FOR EQUATION 2 
X3 : PARTIAL SUM OF EXPECTED NUMBERS IN CELLS FOR EQUATION 3 
X4 : PARTIAL SUM OF EXPECTED NUMBERS IN CELLS FOR EQUATION 4 


GN ER ANGLED ED GY LALIT ID IE IT IE GD IT IE ID ITIL ITD GEIL ID IL IE IT IT GD ITIL ITIL VEN ED IV IVIL IN NGL TL OVAL TOOL GLE GC IN ON ON GL IV ON OS EN GV ON TN GLEN 


TYPE DECLARATION 
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REAL FOM,FOTH,FOV,FOA,P1,P2,P3,P4,Q1,Q2,Q3,Q4, PROD,X1,X2,X3,X4 
RRR RRERERRERERRREREREERERRRRERERERRRREREREERRERRERRERERERRRRRERRERRRER EEE. 
OPENING AN OUTPUT FILE 

CALL EXCMS ('FILEDEF 13 DISK EXPECT DATA Al1') 

SRR ARERRRRRRRERREREICRERRRERRERRRRRERRRERRRRREERERERRRRRERRRRRER 
READING NUMBER OF FAILURES IN EACH TEST INTERACTIVELY 
WRITE(*,*)'PLEASE WRITE NUMBER OF FAILURES IN MANUFACTURER TEST' 
READ(*,*) FOM 

WRITE(*,*) 'PLEASE WRITE NUMBER OF FAILURES IN TEMP.AND HUM.TEST' 
READ(*,*) FOTH 

WRITE(*,*) 'PLEASE WRITE NUMBER OF FAILURES IN VIBRATION TEST' 
READ(*,*) FOV 

WRITE(*,*) ‘PLEASE WRITE NUMBER OF FAILURES IN ALTITUDE TEST’ 
READ(*,*) FOA 
RRERERERERERRERERERREREERRRERIERRERIIIRERERERERRERERERREERRRREREREREERRE 
CALACULATION OF SUCCESS RATIOS 

P1=(20.0-FOM)/8.00 

P2=(20.0-FOTH)/8.00 

P3=(32.0-FOV)/8.00 

P4=(20.0-FOA)/8.00 

RRR RERRERATR ERE ERIEREERKEERERERERERTEEREREREREEREEEEEERRRER 
CALCULATION OF FAILURE RATIOS 

Q1=FOM/8.00 

Q2=FOTH/8.00 

Q3=FOV/8.00 

Q4=FOA/8.00 
PRAIRIE IRE EREEREREEEREERERERREREEERER EE 
CALCULATION PARTIAL SUMS OF EXPECTATION FOR PROGRAM PARAM 

PROD=2 .0*P3+2 .0*Q3+2.0*P4+2.0*Q4 

X1=4.0*P1+4.0*P2+PROD 

X2=4.0*P1+4.0*Q2+PROD 

X3=4 .0*Q1+4.0*P2+PROD 
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X4=4 .0*Q1+4.0*Q2+PROD 
PRR RRR REE RERREEREERREERRERERERRERERRREERERRRREE ER 
WRITING THE RESULTS TO AN OUTPUT FILE 

WRITE (1e,*) X1,X2,X3,X4 

STOP 

END 
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APPENDIX E. PROGRAM PARAM 


RRAERAERERERARERRERRERERERERREERERERERRRRERREBREERRERREREREREREEREERERERE 


PROGRAM PARAM 
HARRKAAKIRERE RARER RERREREREREREREER ERE REE EER EREREREEREERERERR 
THIS IS THE PROGRAM TO CALCULATE PARAMETERS OF LOGLINEAR MODEL 
WITH IMSL SUBROUTINE. IT TAKES PARTIAL SUMS FROM PROGRAM INITIAL 
AND SOLVES FOUR NONLINEAR EQUATIONS, WHICH HAVE FOUR UNKNOWNS. 
THE PROGRAM USES AN IMSL SUBROUTINE CALLED DNEQNF TO SOLVE THIS 
EQUATION.IT WRITES SOLUTIONS OF EQUATIONS TO OUTPUT FILE CALLED 
PARAM DATA. 


REECE SREREEREERRREREERREREREREREEREREERERERERREREREEER 


VARIABLES 


PEREAREREKRAREEREREEREREREREREREEREEREREERERERERERRERERRERERRERRERERERERE 


ITMAX : MAXIMUM ITERATION NUMBER. 

N : PARAMETER 

XGUESS : INITIAL GUESS FOR FOUR NONLINEAR EQUATIONS. 
F : NONLINEAR EQUATIONS. 


KARRI IIIT ARIAT IERIE AIRE RERERREERIREERERREREERI 


TYPE DECLARATION 

PRAIRIE IARI IT AAAI AIA AAA IPRA IARI RIAR RIKER ERR 
PARAMETER (N=4) 

REAL*8 ERRREL 

INTEGER ITMAX,N 


INTEGER K 
REAL*8 FNORM,X(N),XGUESS(N) 
EXTERNAL ACN 


FEENEY AL EL EDEL AE EDEL AL EL GL EL AL AL ED GLAD GLEN ED ISIS IL GL IL GD ILL IN TREN GN IN AU IN EDEN OU IUGR IU OL ON 


OPENING A FILE FOR WRITING RESULTS 
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CALL EXCMS ('FILEDEF 9 DISK PARAM DATA A') 

INITIAL GUESS 

DATA XGUESS/3.5D0,3.5D0,3.5D0,3.5D0/ 

INITIALIZATION 

ERRREL = 0.0001D0 

ITMAX = 10000 

CALLING OF IMSL SUBROUTINE 

CALL DNEQNF ( ACN,ERRREL,N,ITMAX,XGUESS, X ,FNORM) 

GR CIEOT AERO OOS CATS CNL CSTD ICICI CCRT ICICI CORSICA SCI. ICCA TU SCASCA IC ACSI CCCI ICRC OSTA 
RESULTS 

WRITE (9,1) (X(K),K=1,N),FNORM 

FORMAT('THE SOLUTION TO THE SYSTEM IS',/,'X=(',4F8.2,')',/, ‘WITH 
+FNORM=' ,F8.2,//) 

END 


dette te ee eee ee ee eT EE PR eB ER RAH REA ee tee See eset testes 
FREER IL AL GL LAL FL AL GL GY ENGL GL ENGL EDEN GTID GD EL GEED AS GEGVEL TS GL IL ELEY EVAL ED GL GN GLOVE IEE EDEN GL EL EY EV EL IVE EVEL VEL IL ED EYED ED IDED 


SUBROUTINE ACN (X,F,N) 
PHEKAAERARERRERERERERERRRERRERREEERER EEE EERERRERERRER ERE ERE 
VARIABLES - 
X : INITIAL GUESS 
F : NONLINEAR EQUATIONS 
KEKE ERIRERRE ERIE ERE ERR ER ERR ER ERR ERERERRERERERRERRE 
TYPE DECLARATION 
REAL*8 X(N) ,FCN) 
INTEGER N 
RAKE EKIAAKARRAAEEAAEREERRERERREREREREEEKEEREEERERERREREEEREREREERE 
1 ST EQUATION 
F(1)=DEXP(X(1) )*(DEXP(X(2)+X(3)+X(4) )+ 
(1/DEXP(X(2)+X(3)+X(4)))+ 
(DEXP(X(2)+X(3)) /DEXP(X(4)))+ 
(DEXP(X(2)+X(4) ) /DEXP(X(3) ))+(DEXP(X(3) 
+X(4) ) /DEXP(X(2) ))+(DEXP(X(2) )/DEXP(X(3)+X(4)))+ 


+ + + + 
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~ (DEXP(X(3) )/DEXP(X(2) 
+ +X(3)))+(DEXP(X(4) ) /DEXP(X(2)+X(3))))-LOG(32.00) 

BRERA REERERERRRRRERREERRERERRERERERERRERRERRERRRERERE RREREREEEREE 
2 ST EQUATION 

F(2)=DEXP(X(1)+X(2) )*(DEXP(X(3)+X(4) )+ 
+ (1/DEXP(X(3)+X(4) ) )}+(DEXP(X(3))/ 
+ DEXP(X(4) ) )+(DEXP(X(4) ) /DEXP(X(3))))-LOG(23.00) 

Setetetete RRR AEEEREREEERRERERRERRREERRERRRERERRRERERREERERERERRRERRERARRR EERE 
3 ST EQUATION 

F(3)=DEXP(X(1)+X(3) )*(DEXP(X(2)+X(4) )+ 

+ (1/DEXP(X(2)+X(4) ) )}+(DEXP(X(2))/ 

+ DEXP(X(4) ))+(DEXP(X(4) ) /DEXP(X(2))))-LOG(23.00) 

Peete RRR RREEREREREEREREEREREREREEEREREEERRREREIREREREERRERERERREEREREE 
4 ST EQUATION 

F(4)=DEXP(X(1)+X(4) )*(DEXP(X(2)+X(3))+ 
BF (1/DEXP(X(2)+X(3) ))+(DEXP(X(2))/ 
+ DEXP(X(3) ))+(DEXP(X(3) )/DEXP(X(2))))-LOG(14.00) 

RETURN 

END 
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APPENDIX F. PROGRAM LLMDEP 


rales labs aleals aleoleelealsalsaleals oleals als aloals lsolecleals ale alsalsale leaks eloals ale aleaks ale leoteals ale ate aks alools looks ale ateales ale als leolralsals aloalsrelsals lealswls liwleels als 
FR GN GN GN ENGL EVEN GL EVV GL OLIN TD EVEL ED ILD EL IL OL ED SE EL GL ELON EL OV EV ED EVEL IN EDEL ID OIV EVEN FL IVES TNO IL ED UD IL IN ED GE EL EV EVES ES EV ED ON ASOD 


PROGRAM LLMDEP 

THIS IS THE FORTRAN PROGRAM TO CALCULATE THE RELIABILITY 
OF PYROTECHNIC DEVICE. IT ASSUMES THAT THERE IS A DEPENDENCE 
BETWEEN MANUFACTURER AND ENVIRONMENT TESTS. EXPECTATION-MAXIMIZA_ 
TION ALGORITHM IS USED IN THIS MODULE. ALGORITHM STARTS WITH 
INITIAL GUESSES FOR PARAMETERS AND ESTIMATES EXPECTATIONS. IT 
RECALCULATES CELL PROBABILITIES AND UPDATES EXPECTATIONS UNTIL 
IT CONVERGES. AN ITERATIVE NEWTON AND RAPHSON PROCEDURE IS USED 
DURING UPDATATION OF CELL PROBABILITIES. THIS PROCEDURE IS DONE 
BY A SUBROUTINE NAMED UCPROB. 


FRUIT AL AL AEE AT GU AN GU GL GUL GOED INGLY IY GV EY GN GUID IN GU GY ITIL GU GY GE GL ID IN GL IN GD EL GL GL GD ITIVGVINGD IV IRIN GY IRI GL GL GUIDE EV ED EL GD ID 


VARIABLES 


av av avevavay aves Prey EVAL AE AL AL AL AR IR ARAN ER IUGR AU EN ERAN AN GR ANAL AL IN ALIN ERAT AU IN AU GT IL TL GL GU GL ELIE GU GL GEIL ELIT ALIN GY 


FOM : NUMBER OF FAILURES IN MANUFACTURER TEST. 

FOTH : NUMBER OF FAILURES IN TEMP. AND HUMIDITY TEST. 
FOV : NUMBER OF FAILURES IN VIBRATION TEST. 

FOA : NUMBER OF FAILURES IN ALTITUDE TEST. 

MU : OVERALL MEAN. 

LP1 : MEAN EFFECT OF MANUFACTURER TEST. 

EP2 : MEAN EFFECT OF MANUFACTURER TEST. 

LP3 : MEAN EFFECT OF MANUFACTURER TEST. 

LP4 : MEAN EFFECT OF MANUFACTURER TEST. 


TETHA : TWO WAY INTERACTION TERMS 

RHMLE : RELIABILITY OF DEVICE 

MPPPP  : CELL FREQUENCY WITH RESPECT TO TESTS RESULTS. 

MPOPO : SUM OF CELL FREQUENCIES WHICH HAVE PASSED DEVICES 
MANUFACTURER AND VIBRATION TEST. 
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MPOOO : SUM OF CELL FREQUENCIES WHICH HAVE PASSED DEVICES 
MANUFACTURER TEST. 


TP : INITIAL PROBABILITY VECTOR 
FF : UPDATED (FINAL) PROBABILITY VECTOR 
y : CELL EXPECTATION VECTOR 


A,B,C : SOME TERMS TO MAKE THE CALCULATIONS EASY. 
EXPPPP : EXPECTED NUMBER OF DEVICES IN CELL WHICH 
HAS A RESULTANT VECTOR ( P P P P ) IN MANUFACTURER, 
TEMPERATURE, VIBRATION AND ALTITUDE TEST RESPECTIVELY 
FLAG : INDICATOR VARIABLE OF CONVERGENCE FOR PARAMETERS. 
RHMLE : ESTIMATED RELIABILITY OF PYROTECHNIC DEVICE. 


eclecleateclowle tontoate wlowlontealaabowte ake leche abe hohe abe headin ale che wte toaln ale whe ale wlonke alowke ato clealsalealoatvatoats 
FUEL EV GEOR GU EC EC EN DL FU GD ED GN IT IV IN IU UU ILIV EID ID FT ICON EC FD ID TLIC IU IV ID OC IN GL IN EY EDEL OY OD 


ee ee ee era 
CUCU OREN OC IV EN GN ED TN ON ED ED ED EDD 


Feevevouee 


TYPE DECLERATION 

PARAMETER (K=10000) 

LOGICAL FLAG( 16) 

INTEGER I 

REAL FGM; "FOTH, TOVJerOA, LPI, LP2Z 5 EPS, LPS) hhainn 
REAL MPPPP APPPEP Pere seer , 


aE WEPPP APPPEO NEPEP MPR , 
+ MFPPP FPF MePPP Se PEF, 
+ MNPFPP MPP err rae re 
REAL EXPPPP(K) ,EXPPPF(K) ,EXPPFP(K) ,EXPPFF(K) , 
+ EXPFPP(K) ,EXPPPPOK) ,EXPPFP(K) jEXPFFF ORK), 
+ EXFPPP(K) , EXFPPF(K) ,EXFPFP(K) ,EXFPFF(K), 
+ EXFFPP(K) ,EXPPPFC(K) /BXFFFP(K) ,EXPFFFCK) 


REAL MPO000, MNFOOO, MPPOO, MPOPO, MPOOP, "A, B, C 
REAL IPC 6,1); FPC16;2) 920160 

REAL Y(16) 

COMMON / PROB / Y 


vse oe ee UN ON TN ON 


CALL EXCMS ('FILEDEF 15 DISK END DATA Al '' ) 


INITIALIZATION 
RHMLE = 0.0 
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EPS = 0.001 
DO 10 L=1,16 
FLAG(L)=.FALSE. 
CONTINUE 
READING THE NUMBER OF FAILURES IN EACH TEST INTERACTIVELY 
WRITE(*,*)'PLEASE ENTER THE # OF FAILURES IN MANUFACTURER TEST' 
READ(* ,*) FOM 
WRITE(*,*)'PLEASE ENTER THE # OF FAILURES IN TEMP. AND HUM. TEST' 
READ(*,*) FOTH 
WRITE(*,*)'PLEASE ENTER THE # OF FAILURES IN VIBRATION TEST' 
READ(*,*) FOV 
WRITE(*,*)'PLEASE ENTER THE # OF FAILURES IN ALTITUDE TEST’ 
READ(*,*) FOA 
Rated Ra REAR E R RRA ERR ERIE RREER RARER RAE ERRRRRERE 
READING THE INITIAL GUESS FOR EACH CELL IN HYPOTHETICAL 
CONTINGENCY TABLE 
WRITE (*,*)'PLEASE ENTER INITIAL GUESS FOR MPPPP' 
READ(*,*) MPPPP 
WRITE (*,*)'PLEASE ENTER INITIAL GUESS FOR MPPPF' 
READ(*,*) MPPPF 
WRITE (*,*)'PLEASE ENTER INITIAL GUESS FOR MPPFP' 
READ(*,*) MPPFP 
WRITE (*,*)'PLEASE ENTER INITIAL GUESS FOR MPPFF' 
READ(*,*) MPPFF 
WRITE (*,*)'PLEASE ENTER INITIAL GUESS FOR MPFPP' 
READ(*,*) MPFPP 
WRITE (*,*)'PLEASE ENTER INITIAL GUESS FOR MPFPF' 
READ(*,*) MPFPF 
WRITE (*,*)'PLEASE ENTER INITIAL GUESS FOR MPFFP' 
READ(* ,*) MPFFP 
WRITE (*,*)'PLEASE ENTER INITIAL GUESS FOR MPFFF' 
READ(*,*) MPFFF 
WRITE (*,*)'PLEASE ENTER INITIAL GUESS FOR MFPPP' 
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READ(*,*) MFPPP 
WRITE (*,*) ‘PLEASE 
READ(*,*) MFPPF 

WRITE (*,*) ' PLEASE 


READ(*,*) MFPFP 
WRITE (*,*)'PLEASE 
READ(*,*) MFPFF 
WRITE (*,*) ‘PLEASE 
READ(*,*) MFFPP 


WRITE (*,*)'PLEASE 


READ(>*) MPEEF 
WRITE (*,*)'PLEASE 
READ(*,*) MFFFP 
WRITE (*,*)'PLEASE 
READ(*,*) MFFFF 


Petedesbebeshabebetctesbegb oh fetetote micskockosts 
CAI DAT SLR DR RTI E TSI LE KAITAIA ASEAICA SDAA RATS ACR IEEE SCASCRC ACA RACISM AEASLANA SE ASTIS CASE) SUL AT ERTS SERS 


CALCULATION OF INITIAL PROBABALITIES USING CELL FREQUENCIES 


ENTER 


ENTER 


ENTER 


ENTER 


ENTER 


ENTER 


ENTER 


INITIAL 


INITIAL 


INITIAL 


INITIAL 


INITIAL 


INITIAL 


INITIAL 


GUESS FOR 


GUESS FOR 
GUESS FOR 
GUESS FOR 
GUESS FOR 
GUESS FOR 


GUESS FOR 


MFPPF' 


MFPFP' 


MFPFF' 


MFFPP' 


MFFPF' 


MFFFP' 


MFFFF '' 


IP(1,1) = MPPPP/92.00 
IP(2,1) = MPPPF/92.00 
IP(3,1) = MPPFP/92.00 
IP(4,1) = MPPFF/92.00 
IP(5,1) = MPFPP/92.00 
IP(6,1) = MPFPF/92.00 
IP(7,1) = MPFFP/92.00 
IP(8,1) = MPFFF/92.00 
IP(9,1) = MFPPP/92.00 
IP(10,1)= MFPPF/92.00 
IP(11,1)= MFPFP/92.00 
IP(12,1)= MFPFF/92.00 
IP(13,1)= MFFPP/92.00 
IP(14,1)= MFFPF/92.00 
IP(15,1)= MFFFP/92.00 


IP(16,1)= MFFFF/92.00 


ew FEV EVIL OU EN AN GN IL INGE IV ED ED EVEL OV ADIN IV IVIDEV IDV ED GE GT EV IV ED GIVEN GYD ID ID ID ID ED IN ID IV IV ED ED IN GD ID IVIL EV ID IVIV IVIL ITED IV IV INET ON TD OD 


* DETERMINATION OF INITIAL FREQUENCIES FOR LIKELIHOOD ESTIMATION 
Y(1) = MPPPP 
Y¥(2) = MPPPF 
Y(3) = MPPFP 
Y(4) = MPPFF 
Y(5) = MPFPP 
Y(6) = MPFPF 


Y(7) = MPFFP 
Y¥(8) = MPFFF 
C3) = MFPPP 
Y(10)= MFPPF 
Y¥(11)= MFPFP 
Y(12)= MFPFF 
Y¥(13)= MFFPP 
Y(14)= MFFPF 
Y¥(15)= MFFFP 
Y(16)= MFFFF 


# PRAIA IIIA IIR IIIA IRI IIIA ITER IR IRI IRI AIT RI RII ERE 
* CALL A SUBROUTINE WHICH UPDATES CELL PROBABILITIES USING NEWTON 
aS AND RAPHSON PROCEDURE 

CALL UCPROB( IP,FP ) 
* VRPT RAIA IATA AAI API IRI AA AAI AARP EE ERE IRIE IEA AIK IAIREREE 


UPDATATION OF CELL FREQUENCIES 
MPPPP=FP(1,1)*92.00 
MPPPF=FP(2,1)*92.00 
MPPFP=FP(3,1)*92.00 
MPPFF=FP(4,1)*92.00 
MPFPP=FP(5,1)*92.00 
MPFPF=FP(6,1)*92.00 
MPFFP=FP(7,1)*92.00 
MPFFF=FP(8,1)*92.00 
MFPPP=FP(9 ,1)*92.00 
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MFPPF=FP(10,1)*92.00 
MFPFP=FP(11,1)*92.00 

MFPFF=FP(12,1)*92.00 

MFFPP=FP(13,1)*92.00 

MFFPF=FP(14,1)*92.00 

MFFFP=FP(15,1)*92.00 

MFFFF=FP(16,1)*92.00 

Sededetedede ec RRR RERRERRREERERRERERREEEREERERRRERRERERRERRAERERAREEER 
INITIAL EXPECTATIONS 

EXPPPP(1)=MPPPP 

EXPPPF(1)=MPPPF 

EXPPFP( 1)=MPPFP 

EXPPFF( 1)=MPPFF 

EXPFPP(1)=MPFPP 

EXPFPF(1)=MPFPF 

EXPFFP(1)=MPFFP 

EXPFFF(1)=MPFFF 

EXFPPP( 1)=MFPPP 

EXFPPF(1)=MFPPF 

EXFPFP( 1)=MFPFP 

EXFPFF(1)=MFPFF 

EXFFPP(1)=MFFPP 

EXFFPF (1)=MFFPF 

EXFFFP(1)=MFFFP 

EXFFFF(1)=MFFFF 
PORATED EIR EIRRIRERETERERREREREEER 
MP000=MPPPP+MPPPF+MPPFP+MPPFF+MPFPP+MPFPF+MPFFP+MPFFF 
MFO00=MFPPP+MFPPF+MFPFP+MFPFF+MFFPP+MFF PF+MFFFP+MFFFF 
MPPO0=MPPPP+MPPPF+MPPFP+MPPFF 
MPOPO=MPPPP+MPPPF+MPFPP+MPF PF 

MPOOP=MPPPP+MPPFP+MPF PP+MPFFP 

A=92.00-MPP00 

B=92.00-MPOPO 
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C=92 .00-MPOOP 
NEXT EXPECTATIONS 
DO 20 I1=2,K 
EXPPPP( I )=(20.00-FOM)*(MPPPP/MP000)+(20.00-FOTH)*(MPPPP/MPP00 )+ 
(32.00-FOV)*(MPPPP/MPOPO )+(20.00-FOA)*(MPPPP/MPOOP) 
IF(ABS(EXPPPP(1)-EXPPPP(I-1)).LE.EPS) THEN 
FLAG(1)=. TRUE. 
END IF 
EXPPPF(1)=(20.00-FOM)*(MPPPF/MP000)+(20.00-FOTH)*(MPPPF/MPPO0O)+ 
(32.00-FOV)*(MPPPF/MPOPO) +FOA*(MPPPF/C) 
IF(ABS(EXPPPF(1I)-EXPPPF(I-1)).LE.EPS) THEN 
FLAG(2)=.TRUE. 
END IF 
EXPPFP(1)=(20.00-FOM)*(MPPFP/MP000)+(20.00-FOTH)*(MPPFP/MPP0O)+ 
FOV*(MPPFP/B)+(20.00-FOA)*(MPPFP/MPOOP) 
IF(ABS(EXPPFP(I)-EXPPFP(I-1)).LE.EPS) THEN 
FLAG(3)=.TRUE. 
END IF 
EXPPFF(1)=(20.00-FOM)*(MPPFF/MP000)+(20.00-FOTH)*(MPPFF/MPP00)+ 
FOV* (MPPFF/B)+FOA*(MPPFF/C) 
IF(ABS(EXPPFF(I)-EXPPFF(I-1)).LE.EPS) THEN 
FLAG(4)=. TRUE. 
END IF 
EXPFPP(1)=(20.00-FOM)*(MPFPP/MP000)+FOTH*(MPFPP/A)+(32.00-FOV)* 
(MPFPP/MPOPO)+(20.00-FOA)*(MPFPP/MPOOP) 
IF(ABS(EXPFPP(1)-EXPFPP(I-1)).LE.EPS) THEN 
FLAG(5)=.TRUE. 
END IF 
EXPFPF(1I)=(20.00-FOM)*(MPFPF/MP000)+FOTH*(MPFPF/A)+(32.00-FOV)* 
(MPFPF/MPOPO )+FOA*(MPFPF/C) 
IF(ABS(EXPFPF(1)-EXPFPF(I-1)).LE.EPS) THEN 
FLAG(6)=.TRUE. 
END IF 
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EXPFFP(I)=(20.00-FOM)*(MPFFP/MP000 )+FOTH*(MPFFP/A)+FOV*(MPFFP/B)+ 
(20 .00-FOA)*(MPFFP/MPOOP) 
IF(ABS(EXPFFP(I)-EXPFFP(I-1)).LE.EPS) THEN 

FLAG(7)=. TRUE. 
END IF 
EXPFFF(1)=(20.00-FOM)*(MPFFF /MP000)+FOTH*(MPFFF/A)+FOV*(MPFFF/B)+ 
FOA*(MPFFF/C) 
IF(ABS(EXPFFF(I)-EXPFFF(I-1)).LE.EPS) THEN 

FLAG(8)=. TRUE. 
END IF 
EXFPPP( I )=FOM*(MFPPP/MF000 )+FOTH*(MFPPP/A)+FOV* (MFPPP/B)+FOA* 
(MFPPP/C) 
IF(ABS(EXFPPP(1I)-EXFPPP(I-1)).LE.EPS) THEN 

FLAG(9)=.TRUE. 
END IF 
EXFPPF (1)=FOM*(MFPPF/MF000 )+FOTH* (MFPPF/A)+FOV*(MFPPF/B)+FOA* 
(MFPPF/C) 
IF(ABS(EXFPPF(1)-EXFPPF(I-1)).LE.EPS) THEN 

FLAG(10)=.TRUE. 
END IF 
EXFPFP( 1 )=FOM*(MFPFP/MF000)+FOTH*(MFPFP/A)+FOV*(MFPFP/B)+FOA* 
(MFPFP/C) 
IF(ABS(EXFPFP(1)-EXFPFP(I-1)).LE.EPS) THEN 

FLAG(11)=.TRUE. 
END IF 
EXFPFF(1)=FOM*(MFPFF/MF000)+FOTH*(MFPFF/A)+FOV*(MFPFF/B)+FOA* 
(MFPFF/C) 
IF(ABS(EXFPFF(1)-EXFPFF(I-1)).LE.EPS) THEN 

FLAG(12)=. TRUE. 
END IF 
EXFFPP(1)=FOM*(MFFPP/MF000)+FOTH*(MFFPP/A)+FOV*(MFFPP/B)+FOA* 
(MFFPP/C) 
IF(ABS(EXFFPP(1)-EXFFPP(I-1)).LE.EPS) THEN 

FLAG(13)=.TRUE. 
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END IF 
EXFFPF (I )=FOM*(MFFPF/MF000)+FOTH*(MFFPF/A)+FOV*(MFFPF/B)+FOA* 
(MFFPF/C) 
IF(ABS(EXFFPF(1)-EXFFPF(I-1)).LE.EPS) THEN 
FLAG(14)=.TRUE. 
END IF 
EXFFFP( I )=FOM*(MFFFP/MF000)+FOTH*(MFFFP/A)+FOV*(MFFFP/B)+FOA* 
(MFFFP/C) 
IF(ABS(EXFFFP(1)-EXFFFP(I-1)).LE.EPS) THEN 
FLAG(15)=.TRUE. 
END IF 
EXFFFF (I )=FOM*(MFFFF/MF000 )+FOTH*(MFFFF/A)+FOV*(MFFFF/B)+FOA* 
(MFFFF/C) 
IF(ABS(EXFFFF(1)-EXFFFF(I-1)).LE.EPS) THEN 
FLAG(16)=.TRUE. 
END IF 
Rietedetete dete ERERERERERRRRRERIEARARRERERRREREEERREEREREREREREERERERR CI 
MPPPP = EXPPPP(I) 
MPPPF = EXPPPF(I) 
MPPFP = EXPPFP(I) 
MPPFF = EXPPFF(I) 
MPFPP = EXPFPP(I) 
MPFPF = EXPFPF(I) 
MPFFP = EXPFFP(1) 
MPFFF = EXPFFF(1) 
MFPPP = EXFPPP(I) 
MFPPF = EXFPPF(1) 
MFPFP = EXFPFP(1) 
MFPFF = EXFPFF(1) 
MFFPP = EXFFPP(1) 
MFFPF = EXFFPF(I) 
MFFFP = EXFFFP(1) 
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ot 


MFFFF = EXFFFF(I) 
CHECK FOR THE STOPING CONDITION 
IF( FLAG(1).AND.FLAG(2).AND.FLAG(3) .AND.FLAG(4) . AND.FLAG(5) 
+ .AND.FLAG(6).AND.FLAG(7).AND.FLAG(8).AND.FLAG(9) .AND.FLAG( 10) 
+ .AND.FLAG(11).AND.FLAG(12).AND.FLAG(13).AND.FLAG(14).AND. 
+ FLAG(15).AND.FLAG(16)) THEN 
Dt RRR RRRRERREEEERRERERRRRERERARRRRERRRRECERAREEHAAR ECE 
CALCULATION OF FINAL EXPPPP (STOPING CONDITION IS MET) 
RHMLE=EXPPPP(1)/92.00 
GO TO 30 
END IF 


FR GE AE TE GE INGE IVER GETS TE IC EDGY TEIN GU IE GEIL GUID GT IN IN GUILT ID UGE ID ITIL ID GE IVIL ILITINIVIVITIVGICICINGD IDENT ITIL ICI ININININ IN ID GD 


STOPING CONDITION IS NOT MET. PROBABILITIES FOR THE NEXT 
NEWTON AND RAPHSON PROCEDURE 
FEC, 1) MPPPP/92.00 
PPC 2,1) = MPPPr/ 92.06 
TPC3), 1)°="NPPFP/92700 
IP(4,1) = MPPFF/92.00 
TPS, ) = MPEPE/ 92.00 
IP(6,1) = MPFPF/92.00 
TPG? 51) = MPFFP/92-00 
IP(8,1) = MPFFF/92.00 
IP(9,1) = MPPPP/92.00 
IP(10,1)= MFPPF/92.00 
IP(11,1)= MFPFP/92.00 
IP(12,1)= MFPFF/92.00 
IP( 18, = MPFRP/92 .00 
IP(14,1)= MFFPF/92.00 
PP CIS, 1 = tere 9200 
IP(16,1)= MFFFF/92.00 


v camry 


CALL UCPROB(IP,FP) 
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MPPPP=FP(1,1)*92 
MPPPP=PPe2 192. 
MPPFP=FP(3 ,1)*92 
MPPFF=FP(4,1)*92 
MPFPP=FP(5,1)*92. 
MPFPF=FP(6,1)*92 
MPFFP=FP(7,1)*92 
MPFFF=FP(8,1)*92. 
MFPPP=FP(9, 192 
MFPPF=FP(10,1)*92 
MFPFP=FP(11,1)*92 
MFPFF=FP(12,1)*92 
MFFPP=FP(13,1)*92 


MFFPF=FP(14,1)*92. 
MFFFP=FP(15,1)*92. 
MFFFF=FP(16,1)*92. 


.00 


00 


.00 
.00 


00 


.00 
.00 


00 


.00 


.00 
.00 
.00 
.00 
00 
00 
00 


RRAECEREKEEEREREREREREREREREREREREREREEEREREREERERERERRERERERERERREREERE 


MPOOO=MPPPP+MPPPF+MPPFP+MPPFF+MPFPP+MPF PF +MPFFP+MPFFF 
MFOOO=MFPPP+MFPPF+MFPFP+MFPFF+MFFPP+MFFPF+MFFFP+MFFFF 
MPPOO=MPPPP+MPPPF+MPPFP+MPPFF 
MPOPO=MPPPP+MPPPF+MPFPP+MPFPF 
MPOOP=MPPPP+MPPFP+MPFPP+MPFFP 


A=92 .00-MPP00 
B=92.00-MPOPO 
C=92.00-MPOOP 


KREREKEEKERRERERREREREEKRERERERRERREEREREREEREREEREERREEREREEREEREREREREERERE 


CONTINUE 


RREREREREREREEERERERRREREERERERERERERREEREREREEREREEREREREEREREREREER 


WRITE(15,40)FOM,FOTH, FOV, FOA,RHMLE 
RORMAIC) soxen CASE, ,4X 65.2, 2X, F5-2,2%,F5.2,2%,F5.2,/, 15x, 


+ MEE = * FI? .7) 


RRKEARERERRRAEERREREE REIT ICN IS I IST ICIS TOTS TOTS TOUTS TC TCS CTE TOTS TOTTI ICTS STOTT TS PETE ICT TS 


STOP 
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END 


SUBROUTINE UCPROB (IP,FP) 

RANA AAA RAR ARREAE RRR ER ERE RE RR ERE AER ER ER ER RRR EER ERE REREREAER EY 
THIS SUBROUTINE UPDATES CELL PROBABILITIES USING NEWTON AND 

RAPHSON PROCEDURE WHICH IS DESCRIBED IN SAS . 


RRRERERERRERRRERRRRRERRRRRERRERERRERERREERRRRRRREERRRRERRERERRR RE RRRERRRER 


VARIABLES 

deeded ek kk dete RARER REAR ERAN EERE ERR RRR RRR ERED RRR RR 
FOM : NUMBER OF FAILURES IN MANUFACTURER TEST. 

FOTH : NUMBER OF FAILURES IN TEMP. AND HUMIDITY TEST. 
FOV : NUMBER OF FAILURES IN VIBRATION TEST. 

FOA : NUMBER OF FAILURES IN ALTITUDE TEST. 

MU : OVERALL MEAN. 

LP1 : MEAN EFFECT OF MANUFACTURER TEST. 

LP2 : MEAN EFFECT OF MANUFACTURER TEST. 

LP3 : MEAN EFFECT OF MANUFACTURER TEST. 

LP4 : MEAN EFFECT OF MANUFACTURER TEST. 


TETHA : TWO WAY INTERACTION TERMS 

RHMLE : RELAIBILITY OF DEVICE 

MPPPP : CELL FREQUENCY WITH RESPECT TO TESTS RESULTS. 

MPOPO : SUM OF CELL FREQUENCIES WHICH HAVE PASSED DEVICES 
MANUFACTURER AND VIBRATION TEST. 

MPOO0O : SUM OF CELL FREQUENCIES WHICH HAVE PASSED DEVICES 
MANUFACTURER TEST. 


IP : INITIAL PROBABILITY VECTOR 
PE : UPDATED (FINAL) PROBABILITY VECTOR 
ay: : CELL EXPECTATION VECTOR 


A,B,C : SOME TERMS TO MAKE THE CALCULATIONS EASY. 
EXPPPP : EXPECTED NUMBER OF DEVICES IN CELL WHICH 
HAS A RESULTANT VECTOR ( P P P P ) IN MANUFACTURER, 
TEMPERATURE, VIBRATION AND ALTITUDE TEST RESPECTIVELY 
SIGN : INDICATOR VARIABLE OF CONVERGENCE FOR PARAMETERS. 
RHMLE : ESTIMATED RELIABILITY OF PYROTECHNIC DEVICE. 


LZ 


50 


60 
70 


80 


he hee wh ahs whe ahr el whe ake ale ate howls ake whe elie ahs whee ale whe ate ale ale ate atone abe abowle ale atonts ake ate abe ake ale ake leat ate alo ale ate ale atealeatonke ale atowle ale alealeakealeale ate ale wloalealoale 
REN EUAN FR Ae GR AU AE AN ANGE EY IY ER ER EL ENGL ENGLER GU GR EL ELEN IID IN GLAD EL GD ED EY IN GD IL IL EY ED ED IRIN IL IVIV EL EL ELEC EL ELEN EN ELEC ES TS IVIV ESTED 


TYPE DECLERATION 

LOGICAL SIGN(5) 

REAL IP(16,1),FP(16,1) 

REAL FO@IS,1),F1(15,1),S(15,15),x(15,5) ,BOGS,1)¢ 
FPlse ioe 0( 15,1) ,P11C15 91) 96(5,5),6(5,1),SINV(15 984 F( 15,1), 
+PR|(15,1),PR2(5,1),PR3(15,5),PR4(5,5),SUM,LAST,PI01(16,1),LHE, 
+LHEMAX ,RIP(15,1),PR5(15,5),CINV(5,5) ,DIF(15,1),PR6(5,1),LAMBDA, 
+PI11(16,1),EPS,BLAST(5,1),PLAST(15,1),BNEW(5,1),FLAST(15, 1) 
+,U,XT(5,15),PR4INV(5,5) 

REAL Y(16) 

COMMON / PROB / Y 

INTEGER I,J,K 


RERAERAERRAERERRERERARERERERREEREREREREERRERERRRRERREEREREREERREREREREERERE 


INITIALIZATION 


DO 50 I=1,5 

SIGN(I)= .FALSE. 
CONTINUE 
DO 70 I=1,5 

DO'60) J=1,5 

S(I,J)=0.0 

CONTINUE 
CONTINUE 
RikkkehkkeKhkerKeRRR RR Ra LK, RRR RR RR RR RR RR RR RR RR EER 
READING THE DESIGN MATRIX 
CALL EXCMS ('FILEDEF 9 DISK DESIGN INPUT Al") 
DO 80 I=1,15 

REAw@oas) X(1,1),XGl,2),.X(.3) XC 1,09, X(1,5) 
CONTINUE 
REWIND 9 
REREERREREREREEERRARRERRERRERERRRERERERERRERRERERRERREIREREERRRRERERREARRRR RE 
INVERSE OF VARIANCE AND COVERIANCE MATRIX FOR INITIAL BO 
DOmO0 tts 


Ue 


100 


DO 90 J=1,15 
IF(1.EQ.J) THEN 
SINV(I,1)=(IP(I,1)-(IP(I,1)#*2.0))*92.00 
END IF 
IF(I.NE.J) THEN 
SINV(I,J)=(-IP(J,1))*IP(1,1)*92.00 
END IF 
CONTINUE 
CONTINUE 
sededetetectedotededs dete eR ERASER RREREREREREERERAREEREREEESR 
LOGIT RESPONSE FUNCTIONS 
DO 110 I=1,15 
F(I,1) =ALOG(IP(I,1)/IP(16,1)) 
CONTINUE 
RAKE RRR EIRENE ERR EEE 
THE TRANSPOZE OF THE DESIGN MATRIX 
CALL TRNRR (15,5,X,15,5,15,XT,5) 


PERRIER 


MATRIX MULTIPLICATION PR1=(SINV*F) 

CALL MRRRR (15,15,SINV,15,15,1,F,15,15,1,PR1, 15) 

PRIA EREARERETAAERIARRERIEEEERRETAEREREEEERERERERRERRREE 
MATRIX MULTIPLICATION PR2=(XT*PR1) 

CALL MRRRR (5,15,XT,5,15,1,PR1,15,5,1,PR2,5) 


RK KKK KE 


MATRIX MULTIPLICATION PR3=(SINV*X) 

CALL MRRRR (15,15,SINV,15,15,5,X,15,15,5,PR3,15) 

PRR RERKREER ERITREA AER ASIARA EERIE RE 
MATRIX MULTIPLICATION PR4=(XT*PR3) 

CALL MRRRR (5,15,XT,5,15,5,PR3,15,5,5,PR4,5) 

BPAAAASACAS TCR TACHA CARER AASI ARIES SAREE ERIE RRR RRERES 
INVERSE OF THE MATRIX MULTIPLICATION PR4=PR4INV 

CALL LINRG (5,PR4,5,PR4INV,5) 


tiie wae 
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120 


130 


140 


150 


INITIAL ESTIMATION OF PARAMETERS BO 
MATRIX MULTIPLICATION BO=(PR4INV*PR2) 
CALL MRRRR (5,5,PR4INV,5,5,1,PR2,5,5,1,B0,5) 
GOW Ge TCE PETC I IS TOTO TCS TCT TS TCI TCICIS TOTS TOE WINTON TS TO TCTC EC TRICE TRIS TR TS TUNIC TUTE RRRERERER RRR: 
FO= X * BO 
CALL MRRRR (15,5,X,15,5,1,B0,5,15,1,F0,15) 
CASS RSS RATATAT STIR CITI De CTA Be ISCAS LSTA STRUTT ES CIO RETIRE COTES CCRC CCAS TERT CC RET CTCL TL MCSE CSCC 
INITIAL PROBABILITIES PI0=(EXP(FO) ) 
DO.120' 1=1, 15 
PIO(I,1)=EXP(FO(1I,1)) 
CONTINUE 
SUM=0 . 0 
DOws0et=1,15 
SUM=SUM+PI0(I, 1) 

CONTINUE 
LAST=1.0/(1.0+SUM) 
PRAIRIE EERIE ERE EREREREEEREREREREREREERERERERETEERER 
PROBABILITY MATRIX WHICH INCLUDES 16 VALUES PIO1 
FOR THE INITIAL ESTIMATE OF LIKELIHOOD ESTIMATION 
DO 140 I=1,15 

PIO1(1,1)=PI0(I,1)*LAST 
CONTINUE 
PI01(16,1)=LAST 
PREECE RREEREREEEREREREERERREEEERERERERERERERERREREREREREREREREEREREER 
INITIAL LIKELIHOOD FOR NEXT ITERATION AT STEP BO 
LHE=0 . 0 
DO 150 I=1,16 

LHE=LHE+Y(1I)*ALOG(PIO1(I,1)) 
CONTINUE 
LHEMAX=LHE 
RAAT RARER ERCRERERERICRREERREREEERRERREERERREREERERREEREERERERERREREE 
REORGANIZED INITIAL PROBABILITIES FOR UPDATATION RIP 
DO 160 I=1,15 

RIP(I,1)=IP(I, 1) 
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160 


190 


CONTINUE 
CASTOR SIL SCI TSI SCTSERIA SLICES CO CLOSES RSC RICANS DECASCIS CISTI CAISSON TIS CIS CISTI IS CISC CSE TIC CS OS Co CT 
FIRST ITERATION IN NEWTON AND RAPHSON METHOD 
DO 170 l=ieas 
DO 180 J=1,15 
IF(I.EQ.J) THEN 
SINV(I,1)=(PIO(I,1)-(PIOCI,1)#*2.0))*92.00 
END IF 
IF(I.NE.J) THEN 
SINV(I,J)=(-PIO(J,1))*PIO(I,1)*92.00 


END?) TE 
CONTINUE 
CONTINUE 


lhe afhee ahe whe, Soothe ate 
PRE RREERREREEEREREERERERRERREE RREREEREERERRERERREREERERERRERERERE 


MATRIX MULTIPLICATION PR5=(SINV*X) 
CALL MRRRR (15,15 .SINV j15., 15.5.4, 15, 15,5 ,PR55190 


le legocloctecloctoatoate Sefetedodevedcedkede tefasteste leslentedlonkeclecloate. oy eieatactoctanswetoctaesoetseknstat 
PIAA AIEEE RAEERERERREEEEREIEREREEREEREEREREEEREREEEREREEREREERE 


MATRIX MULTIPLICATION C=(XT*PRS5) 
CALL MRRRR (5,15,XT,5,15,5,PR5,15,5,5,C,5) 
PEPTIDASES IKERE EEE RREREREREREEREREREREER 
INVERSE OF THE MATRIX C=CINV 
CALL LINRG (5,C,5,CINV,5) 
PEPRETCH CITC AAA HAC TCA AAC AH PPT ACHAEA ACHAEA AEE AA RETR EERIE EHR ERE 
DO 190 I=1,15 

DIF(I,1)= 92.00*(RIP(I,1)-PIO(I,1)) 


CONTINUE 
PHAKIC AHR AIEEE AAAI RIA ER EERIE RE ER ERERERERE 


MATRIX MULTIPLICATION FOR G=(XT*DIF) 
CALECURRER (5, 15,X1 ,5, [53a WIE, £540, 1,6.5) 


Titictdctcsted 


MATRIX MULTIPLICATION (PR6= CINV*G) 
CALL MRRRR (5,5,CINV,5,5,1,G,5,5,1,PR6,5 ) 


HAKKAR ERERAERERERERIRE REAPER RIKI ERIE ER ERE EERE 
LAMBDA=1.0 
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200 


220 


230 


240 


250 


DO 210 I=1,5 
PR6(I,1)=PR6(1,1)*LAMBDA 
CONTINUE 


whe ake ake afeake ale ahs ake ake aleatenkeake ate ake ale aloats atamls ake nlsaks heals ale ale afeale ate nke ale akeake afealeafeale ate aks ale nle ake ake aks ale ale aleats ale alewks ahealealeale ale ahealealeate ale wa ale 
JEON EN EN PUAN UN AN AU EN ENN EV GN AN GN AU EE GUILE GN IY IN INGE IU GLEE EDEN IN IN GD ANGE IU IU IN GU GE TEIN GV ELEN OL EV IN GL IV GL ENED INI EDEL GD EN OLIN OS 


INITIAL VALUE FOR Bl 
DO 220"1=1,5 

B1(1I,1)=BO(I,1)-PR6(I,1) 
CONTINUE 
MATRIX MULTIPLICATION F1=(X*B1) 
CALL MRRRR (15,5,X,15,5,1,B1,5,15,1,F1, 15) 
CALCULATION OF PROBABILITIES FOR Bl PI1=EFP(F1) 
DO 230 I=1,15 

PI1(1,1)=EXP(F1(1,1)) 
CONTINUE 
PPR RREEREREREREEEEREERERRRERRERERRERERRRREERERRRERREREEREER 
CALCULATION OF THE 16 TH PROBABILITY VALUE 
SUM=0.0 
DO 240 I=1,15 
SUM=SUM+PI1(I, 1) 

CONTINUE 
LAST=1.0/(1.0+SUM) 
SIO CRS CRS A CRIS SLITS CAIUS OY CESSES ORCAS CISCO CG RIERA ACAI AIRISASA LIC RA RCA EAI RA ITCRA CI CO CAA 
DO 250 I=1,15 

PI11(1,1)=PI1(I,1)*LAST 
CONTINUE 
PI11(16,1)=LAST 
SOOO TISCO COSTCO TREO OO USCIS OCS ICICI DCCA CIRO IAI SCAST ANCL LIST RICA SEIS IELCICA PRR 
INITIAL LIKELIHOOD ESTIMATION 
LHE=0 .0 
DO 260 I=1,16 

LHE=LHE+Y(I)*ALOG(PI11(1,1)) 
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260 CONTINUE 


whechretealealeateatsabeals abe aleatoats abe whe whes whe ate whe whe whe whe ale beaks abe ake abe nlowtes whe whe ake ateatoateleatenle ale alewte wheats ale atoatonte he whe ahs whoukc mle whe alowte wloakes ate abealsate ate 
FEIT AT ED EN TL GC ID TD GD DIL GT IY ITED ID EC AN FEV GL AVIV IDET IV FD FD FY IT ID IDE ERED FED ID EY ID ITAL EVID GV IS OD IL IU GV EVID IVIV EV AVES TL EDEN EYED OD 


K=0 
EPS=0.001 
IF(LHE.LT.LHEMAX) THEN 
K=K+1 
IF(K.GT.10) THEN 
DO 270 I=1)5 
BLAST(I,1)=BO(I,1) 
270 CONTINUE 
DO 280 I=1,15 
PLAST(I,1)=PI01(I,1) 
280 CONTINUE 
GO TO 490 
END IF 
LAMBDA=LAMBDA/2 .0 
GO TO 200 
END IF 
LHEMAX=LHE 
DPPC RRR IKEA RRR ERIE EERE EREEREREREREREERRERRRERER 
DO 290 I=1;5 
IF(ABS(BO(I,1)-B1(I,1)).LE.EPS) THEN 
SIGN(I)=.TRUE. 
END IF 
290 CONTINUE 
Ge ee oe EAL IEA RRI TEAST LIAM EAR AIA SCA ICTR OTA ICA ITAR AC EHCIE ORLA ICA FCASTAOT ALA TEATE ISH ASK KIC ESTA MEAITANT AS I AICA IEE SLES AOC RICERCA 


CHECKING CRITERIAS 
IF(SIGN(1).AND.SIGN(2).AND.SIGN(3).AND.SIGN(4).AND.SIGN(5)) THEN 
DO 300 I=1,5 
BLAST(I,1)=B1(I,1) 
300 CONTINUE 
DO 310 I=1,15 
PLAST(I, 1)=P111@0; 1) 
310 CONTINUE 
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320 


+ 


+ 


GO TO 490 
END IF 
Pe Or SOO OK OO CORI URCOR ICO ROCIO CACO DCO OIC ICR DCO OOK IC ORCC DRACO CARA 
CRITERIAS ARE NOT MET THEN NEW ITERATIONS 
pO 320 T=4,5 
BLAST(I,1)=B1(I,1) 
CONTINUE 
DO 330° 1=1/,15 
PLAST(I,1)=PI1(I, 1) 
CONTINUE 
RARER EERE ERR RERRERREREERRERREEER RRRERREEREERERERREEREEERE 
DO 360 I=1,15 
DO 350 J=1,15 
IF(1.EQ.J) THEN 
SINV(I,1)=(PLAST(I,1)-(PLAST(I,1)**2.0)) *92.00 
END IF 
IF(I.NE.J) THEN 
SINV(1,J)=(-PLAST(J,1))*PLAST(1,1)*92.00 
END IF 
CONTINUE 
CONTINUE 


SUES FTES EL GS IU aN aN aU AN TEES OL AU UGE IN GD IE IT ING 


ole. wlowtoale toate 
RERERERRERRRRRERERERREREREERRERERERERRER 


MATRIX MULTIPLICATION PR3=(SINV*X) 

CALL MRRRR (15,15,SINV,15,15,5,X,15,15,5,PR3,15) 
dedetetededeteReteTeteHREEKEERRRRIEREEERRREREERKRERRREERERREREEERRERRERREEREEREEEER 
MATRIX MULTIPLICATION C=(XT*PR3) 

CALL MRRRR (5,15,XT,5,15,5,PR3,15,5,5,C,5) 
REREKEREEERERREEEEEEREREREEKEEERERRERERERRERERERERRRRRERREEERREREREREE 


CALL LINRG ( 5,€,5,CINV,S ) 


oaks lelealeoletroleatsateatoaleutowte. hele boats lealrals beaks at. heakeloalealoal. lealrale ale ate ole ole oaks cloaks ule leaks look leale bl Le oak ob abe eat 
FR ER AEE IE AR AN AR EN EN AR AN AE AN EN GN AN AN EN INGE AN EN AE AN GE AN AE ER AL EL EE AE AN AU Ge aU AE GN EN EN IE TUES ES TUTE GE IS IN IE ID IE IY ID GN GD LIL IE IN YY SD 


DO 370 1[=1,15 
DIF(I,1)= 92.00*(RIP(I,1)-PLAST(I,1)) 
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370 


400 


410 


420 


CONTINUE 

MATRIX MULTIPLICATION G=(XT*DIF) 

CALL MRRRR (5,15,X1,5, 1554, Dik | 1555715055) 
MATRIX MULTIPLICATION PR6=(CINV*G) 

CALL MRRRR (5,5,00NW>5,5,1,6,5,5, 1,286.5) 


he cheeks hoeks elrelrelrale clrelreleatrats cleuleals aloals aboot ate alraleals alealents oleate alealsatoale ale als elealealeals lealealealsalealoat ale alralelealeale leelealealeel lel eleeleeleels 
FR GNONEN GN EN GLEN EC GL ELIE SE TE ID EN ODED ED TD ORGS TEEN ED ED DID ED TU GD IE GE ED IE TCG ICID IDES IL ICED GUID ISIN ID ICID IN ED IDEN EN IN EC ESOL ENED ES OD 


LAMBDA=1.0 
NEW PARAMETER ESTIMATES 
DO 390 I=1,5 

BNEW(I, 1)=BLAST(I, 1)-(LAMBDA*PR6(I,1)) 
CONTINUE 
tevedetotctcdededetdotetedete de RRR RRREREEERERRRERERERRRERRREERARR RARER RE 
MATRIX MULTIPLICATION F1=(X*BNEW) 
CALL MRRRR (15,5,X,15,5,1,BNEW,5,15,1,F1,15) 
DO 400 I=1,15 

PI1(I,1)=EXP(F1(I,1)) 
CONTINUE 
Sosetete teste de ote se eee ee ie esoteric tarwet ee ecu we ae oe 
SUM=0.0 
DO 410 I=1,15 
SUM=SUM+PI1(I,1) 

CONTINUE 
LAST=1.0/(1.0+SUM) 
PRTC CITC TCR AAICRI TRI IIRIRRIREEERREICRERERERRRERERRRRE EEE 
DO 420 I=1,15 

PI11(1,1)=PI1(I,1)*LAST 
CONTINUE 
PI11(16,1)=LAST 
BRR RRR BERR RRER EERE ERERERRERRERERE EERE REE 
LIKELIHOOD ESTIMATION 
LHE=0 .0 
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430 


450 


460 


DO 430 I=1,16 
LHE=LHE+Y(1)*ALOG(P111(I,1)) 
CONTINUE 
RTE RRR RERRRERRRRRRERREERERRERERARRERERRRRRERERERRR REE 
K=0 
EPS=0.001 
IF(LHE.LT.LHEMAX) THEN 
K=K+1 
IF(K.GT.10) THEN 
GO TO 490 
END IF 
LAMBDA=LAMBDA/2.0 
GO TO 380 
END IF 
LHEMAX =LHE 
DO 440 I=1,5 
IF(ABS(B1(I,1)-BLAST(I,1)).LE.EPS) THEN 
SIGN(1)=. TRUE. 
END IF 
CONTINUE 
thkkhikkkeeeeeks 
IF(SIGN(1).AND.SIGN(2) .AND.SIGN(3) .AND.SIGN(4).AND.SIGN(5)) THEN 
DO 450 I=1,5 
BLAST(I,1)=BNEW(I, 1) 
CONTINUE 
DO 460 I=1,15 
PEASTOI 1) =PI11(1, 1} 
CONTINUE 
GO TO 490 
END IF 


FEIN GL CUAL GY IRAN AV EN IVAN IL GL GL IVIL OL GY IY IY IY ID IV EV ED IY IV GY IVEY IETS GV IV EV IC ETAL IV ITIL IV IV GT IVIV ISIC IV IVIL ITIVE INTL IVES TID IN IN OD 


ephcode orbs oles he. 
LIRA RERRREERERRRERERREREERRERERRERERERRR 


CRITERIAS ARE NOT MET THEN NEW ITERATIONS 
DO 470 I=1,5 
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470 


480 


490 


500 


510 


520 


BLAST(I,1)=Bicl a) 
CONTINUE 
DO 480 I=1,15 
PLAST( 151) =Pllucie 
CONTINUE 
GO TO 340 


BERRA EKRERRREREREREREREREREERERERERRERERERERERERERERRERERROCERRE 


CONTINUE 
MATRIX MULTIPLICATION FLAST=(X*BLAST) 
CALL@MERRR (15 ,5,%, 15,51, 8UAST 5,15, 2, FuASt iS) 
HAA EA ERA AAA AAA RRA A ARR RE ERA RARER REBAR RRA RRR EEE 
CALCULATION OF FINAL PROBABILITIES 
DO 500 I=1,15 

PLAST(1I,1)=EXP(FLAST(I,1)) 
CONTINUE 
KAKKERKER KIER AERERERKERIER ARERR ERERERERER REERE BERR RRR REERERERIR 
SUM=0.0 
DO S1051=1,15 

SUM=SUM+PLAST(I,1) 

CONTINUE 


REE SEERA ERECREEEEREEEERERERRREREEREEERE 


LAST=1.0/(1.0+SUM) 

DO) 520. 1=115 
FP(I,1)=PLAST(I,1)*LAST 

CONTINUE 

FP(16,1)=LAST 

RETURN 

END 


lZ2 


APPENDIX G. PROGRAM MLEB 


eho ale ale ols alools ale alsalsaloats als aloals alealsalraloalools ole nls heals aloalsclontcoloaloalsealealoalouts les aleolouls aloolools ale alealools aloals ols elrelsals alealoelosls leelselseleelbaeloeloals 
FEIN ANON ENTS POISED OS OU IN GOOD ON ODED OE ODES TU ETI EN ER ODOR FEIN INGEN GSES IN GL ED ED IDEN FV ICIS ED OLED OD OS ER LON FLED ERED PCED EN ER EN OR ED ON ES 


PROGRAM MLEB 


howls aloals alealsntsaloalsalouts atrale aloalsaloaloalsalealoaloalouls aloalealralonlsaloule alealsulsaloalouloatoals als ale lealoatsaloaloulonlaalsalealoalouls aloaloaloalsalealsaloolontentsaleals 
FUP ON ON TEED IE IL EN INGE GUID EN FHI FEIN IE ID INTUTE IN IY ON GOED FEIN GD ID DINED FD EDIE IE ID IMIS IE IT IN ID ENTE TE IEOD OD IN INTL IN IDEN IN FOILED EY OE 


THIS IS THE PROGRAM TO CALCULATE RELIABILITY OF THE DEVICE 


WITH DEPENDENCE ASSUMPTION.PROGRAM ASSUMES THAT FAILED ITEM FROM 
ANY OF ENVIRONMENT TESTS FAILS FROM MANUFACTURER TEST TOO. THIS 
IS WORST CASE SCENARIO. IT READS NUMBER OF SUCCESSFUL ITEMS FROM 


AN INPUT DATA CALLED SUSVECT. FINALLY THE PROGRAM WRITES RESULTS 


TO AN 


OUTPUT FILE CALLED RESULT. 


whens aloalealoalesalsaloaloulsaloalealealealvalealoaloalesaleatea nla nlewles ale ale ata ale oloulonlas lowloule ale loatsaloaloulonloults alouleuloaloatoalts altoulontcale ale aloulealsalealoalsaloulesaloultsctouts 
FETC OV ON PC OVEN OV EDT ED TS ONIN GOOD ERTL OD IU GE IN ITIL TIN INES TEIN OD IN OD GL GOIN EOIN ONIN OD IL EDEL OE ED IO IN IS ON UU GN TY OD IN ID IN FETE TE ID ID OU IL OD TE 


VARIABLES 


PEPTIC IDIOTS OER EL TEIN DPC IEEE INI IN ID IN GEIL GD ED IO IN IN GD IN GE IE GT ES TT GD IID IVIL IV GV OD IT ICIC IVIL INTL DID ID ITD IN IN IN IV ED OD IV EDEN OD 


SOM 
SOTH 
SOV 
SOA 
R1H 
R2H 


: NUMBER OF SUCCESSFUL ITEMS IN MANUFACTURER TEST 

: NUMBER OF SUCCESSFUL ITEMS IN TEMP. AND HUMIDITY TEST 

: NUMBER OF SUCCESSFUL ITEMS IN VIBRATION TEST 

: NUMBER OF SUCCESSFUL ITEMS IN ALTITUDE TEST 

: ESTIMATED PROBABILITY OF PASSING FROM MANUFACTURER TEST 
: ESTIMATED PROBABILITY OF PASSING FROM TEMPRATURE AND 


HUMIDITY TEST 


: ESTIMATED PROBABILITY OF PASSING FROM VIBRATION TEST 
: ESTIMATED PROBABILITY OF PASSING FROM ALTITUDE TEST 
: ESTIMATED RELIABILITY OF ITEM AFTER EXPOSURE TO 


SEVERAL ENVIRONMENT TESTS. 


: DUMMY VARIABLE 


BRA AERER EERE AR RENE REENE BRRERRCERREREE ERIE EREREREEERERRERRRE 


TYPE DECLARATION 
REAL SOM(5000) ,SOTH(5000) ,SOV(5000) ,SOA(5000) ,X(4) ,R1H,R2H,R3H, 
+ R4H,RHMLE(5000) 
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INTEGER I 
FILES FOR READING AND WRITING 
CALL EXCMS ('FILEDEF 9 DISK SUSVECT DATA A1') 
CALL EXCMS ('FILEDEF 17 DISK RESULT DATA A1') 
READING NUMBER OF SUCCESS IN EACH TEST 
DO 10 I=1,5000 
READ(9,*) SOM(I),SOTH(I),SOV(I) ,SOA(I) 
Seteteteke ese RRRREERRRRREERRERRRRERRERRRRERRRRRRRRERRRERRRRRRRRRRR 
NUMBER OF SUCCESSES IN EACH TEST 
X(1)= SOM(1I)+SOTH(1I)+SOV(1)+SOA(I) 
X(2)= SOTH@) 
X(3)= SOV(I) 
X(4)= SOA(I) 
Dedede detect eR RRR ERRERERERERERERERREREEEREREREREREERERRE 
CALCULATIONS SUCCESS PROBABILITIES IN EACH TEST 
R1H= X(1)/184.0 
R2H= X(2)/40.0 
R3H= X(3)/64.0 
R4H= X(4)/40.0 
RHMLE(I)= R1H * R2H * R3H * R4H 
FERIA AAA ARR AAA RTARATAERERIREIR EREECKEREEERERERRERERERERERERE 
WRITING RESULTS TO AN OUTPUT FILE CALLED RESULT 
WRITE (17,*) RHMLE(I) 
CONTINUE 
STOP 
END 
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APPENDIX H. PROGRAM BONUS 


whe whe ake aheals abe ake ols ols ole ole abs olke ole ale aks ale cleats ole abs ale als ale aleals cleols ols ule ules als ale uke uls ale aks akin) cols ules ale als ale ule ule ale als abe ale leaks ale ole ole als cloaks heels eleuleals cle 
TEGN IVAN SEEN TEIN ED GOITER TE IE GD GO GE TET TL IU TL ED EC IN IR AV EV IV ED ED SD EVID GD ID GD ID ED GD ID IU FEED ID IL IDI ID TE ID GU TDED IN TV OC EN GLOVES GLEN ES 


PROGRAM BONUS 

THIS IS THE PROGRAM TO CALCULATE BONUS PERCENTAGE OF ANY FIRM 
WHOSE LONG RUN SUCCESS PROBABILITIES ARE KNOWN. IN THIS PROGRAM 
IT USES 2000 SUCCESS VECTORS, WHICH ARE GENERATED BY RANVEC IN 
APPENDIX B. THEY ARE GENERATED BY KNOWN LONG RUN PROBABILITIES 
THE PROGRAM USES TWO DATA SETS.THEY ARE PRECALCULATED LCB'S SETS 
-FIRST DATA REPRESENTS FIRST INSPECTION, SECOND DATA REPRESENTS 
SECOND INSPECTION LOWER CONFIDENCE BOUNDS. 

SUCCESS VECTORS REPRESENT OFFERED LOTS.IT HAS A DETERMINISTIC 
BONUS LINE.PROGRAM CALCULATES LCB OF OFFERED LOT,WITH FIRST DATA 
AND COMPARES IT WITH LCB OF BONUS LINE.IF FIRM LCB IS GRATER THAN 
FIRM GETS BONUS. OTHERWISE FIRM HAS A CHANCE TO ONE MORE TRY. IN 
SECOND TRY, PROGRAM CUMULATES SUCCESS VECTORS AND IT USES SECOND 
DATA TO FIND OUT LCB OF CUMULATED LOT. AFTER THIS CALCULATION IT 
COMPARES AGAIN. FINALLY IT COUNTS NUMBER OF TIMES THAT THE FIRM 
GETS THE BONUS IN 1000 REPLICATIONS AND ESTIMATES BONUS PERCENT. 
IT WRITES RESULTS TO AN OUTPUT FILE CALLED BONUS DATA. 


PPAR REACH REAR RERREREREERRERERERERERRERRRRERERREREREREREREREE 


VARIABLES 


BRRERRRRRRRREECREREREREREREREEEREERRERERERRERERERERERERERRERERREREE 


SOM : NUMBER OF SUCCESS IN MANUFACTURER TEST. 

SOTH : NUMBER OF SUCCESS IN TEMPERATURE AND HUMIDITY TEST. 
SOV : NUMBER OF SUCCESS IN VIBRATION TEST. 

SOA : NUMBER OF SUCCESS IN ALTITUDE TEST. 

BLINE : LOWER CONFIDENCE BOUND OF BONUS LINE 

LCB : LOWER CONFIDENCE BOUND 

FILCB : LCB VALUES ARRAY IN FIRST INSPECTION 
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SELCB : LCB VALUES ARRAY IN SECOND INSPECTION 
A,B,C,D : DIMENSIONS FOR USE OF LCB DATAS 


BFI : NUMBER OF TIMES THAT FIRM GAT BONUS AFTER 1 ST INSP. 
BSI : NUMBER OF TIMES THAT FIRM GAT BONUS AFTER 2 ST INSP. 
BTOT : TOTAL NUMBER OF TIMES THAT FIRM GAT BONUS. 

PRCT : BONUS PERCENT. 

COUNT : COUNTER FOR 1000 REPLICATIONS. 

SIGN : INDICATOR OF ACCEPTANCE FOR FIRST INSPECTION 

FLAG : INDICATOR OF ACCEPTANCE FOR SECOND INSPECTION 


on on OS ek a ales whee akrakis ake akes ahs ales wheal afew alee afer abe als alent als aleats on Oe ee a A 
BRERA RRRERRERRRRRRREER RARE RARER RRR ER ERECT BNA ERR CR ERR 


TYPE DECLARATION 
LOGICAL SIGN(4) ,FLAG(4) 

REAL SOM( 2000) ,SOTH(2000) ,SOV(2000) ,SOA( 2000) ,BLINE,LCB,BFI,BSI, 
+ BIOT,FILCB(2,2,3,2) ,SELCB(3,355,4),A,550.0.PRGT COUNT 

INTEGER 1,J,K,L 

PeteTe RRR RRERERERERREEREERERERERERRERERERACERERRREREEREREERERRERRERE 
OPENING FILES FOR READING AND WRITING 

CALL EXCMS ( ' FILEDEF 7 DISK SUCVECT DATA Al‘ ) 

CALL EXCMS ( ' FILEDEF 8 DISK FIRST DATA Al' ) 

CALL EXCMS ( ' FILEDEF 9 DISK SECOND DATA Al‘ ) 

CALL EXCMS’ ( ' FILEDEF 15 DISK BONUS DATA Al' ) 


EV IVIL EV CD IR aN AE AN ee ae ON ae PIV EV AL EVEL IL ERIE IL IL AL IL GS IL EL GD IL IL ID AL EL AL UD EL IL AL IL IN EL IN GD IL IL ER IL IN IN AL IN CLAN IN EN IE GR ER AN AN 


INITIALIZATION 
COUNT = 1.0 
BFI = 0.0 
BSI = 0.0 
BLINE = 0.9250000 
DO 10 I=1,4 
SIGN(I)=.TRUE. 
FLAG(I)=.TRUE. 
CONTINUE 
RRA RARER ARERR ERE REE RERREREREERERREERERERERER IR 
READING FIRST INSPECTION LOWER CONFIDENCE BOUNDS FROM DATA FILE 
DO 50 I=1,2 
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20 
30 
40 
50 


60 
70 
80 
90 


DO 40 J=1,2 
DO 30 K=1,3 
DOr 20" L=12 
READ(8 ,*)eFILCB(1,J,K,L) 
CONTINUE 
CONTINUE 
CONTINUE 
CONTINUE 
READING SECOND INSPECTION LOWER CONFIDENCE BOUNDS FROM DATA FILE 
DO 90 I=1,3 
DO 80 J=1,3 
DO 70 K=1,5 
DO 60 L=1,3 
READ(9,*) SELCB(I,J,K,L) 
CONTINUE 
CONTINUE 
CONTINUE 
CONTINUE 
KEERRERERERIEEREERERREREREREREEEREREEEREEERERRERERERERERERERERREERRE 
READING SUCCESS PROBABILITY OF FIRM IN EACH TEST 
WRITE(*,*) ‘WRITE THE PROBABILITY OF SUCCESS IN MANUFACTURER TEST' 
READ (*,*) PSIM 
WRITE(*,*) 'WRITE THE PROBABILITY OF SUCCESS IN TEMPERATURE AND 
+HUMIDITY TEST' 
READ (*,*) PSITH 
WRITE(*,*) 'WRITE THE PROBABILITY OF SUCCESS IN VIBRATION TEST' 
READ (*,*) PSIV 
WRITE(*,*) 'WRITE THE PROBABILITY OF SUCCESS IN ALTITUDE TEST' 
READ (*,*) PSIA 
REPRE REECREREREERERERREREEREERREREREEERERRERERERERRERRREERERERERRREREER 
READING NUMBER OF SUCCESFUL ITEMS 
FOR FIRST INSPECTION AND SECOND INSPECTION 
DO 100 I=1,2000 


127 


READ(7,*) SOM(1I) ,SOTH(I) ,SOV(I),SOA(I) 
100 CONTINUE 


fe REAR RIAURR ERA EERERAR IRR ERR ERE EERE RRR EERE REER EERE REAR 
sh INSPECTIONS BEGIN 

DOw120 N=1,2000,2 
= HARARE ARE ER RR REE RRR RE RERERRAR RRR ERR RR RARER 
a SUBSCRIPT DEFINITION FOR FIRST INSPECTION 


A=INT(20.0-SOM(N))+1 
B=INT(20.0-SOTH(N))+1 
C=INT(32.0-SOV(N))+1 
D=INT(20.0-SOA(N))+1 


fe Pedder terete RRR EERE ER EREEE RERERERERRERERERRRR ERE RRR 
is CHECK FOR ACCEPTANCE OF FIRST OFFERED LOT 


IF (A.GT.2)) THEN 
SIGN(1)=.FALSE. 

END IF 

Te) (CB .-Cle2)) TREN 
SIGN(2)=.FALSE. 

END IF 

IF (C.GT.3) THEN 
SIGN(3)=.FALSE. 

END IF 

TF (D-GT.2) THEN 
SIGN(4)=.FALSE. 

END IF 

IF(.NOT.(SIGN(1).AND.SIGN(2).AND.SIGN(3).AND.SIGN(4))) THEN 
SIGN(1)=. TRUE. 
SIGN(2)=.TRUE. 
SIGN(3)=.TRUE. 
SIGN(4)=.TRUE. 


GO TO 112 
END IF 

cae REREERERREREEREREREREREEEERRETRECERRCRR CRI IRR RRERRER ERR RRR EER RRR 

* AFTER FIRST INSPECTION DETERMINATION OF LCB OF FIRM 
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LCB=FILCB(A,B,C,D) 
IS LOT WORTH WHILE FOR GETTING BONUS ? 
IF ( LCB.GT.BLINE ) THEN 
BFI=BFI+1.0 
GOmO 110 
END IF 
SUBSCRIPT DETERMINATION FOR SECOND INSPECTION 
A=A+INT(20.0-SOM(N+1)) 
B=B+INT(20.0-SOTH(N+1)) 
C=C+INT(32.0-SOV(Nt1)) 
D=D+INT(20.0-SOA(Nt1)) 
PRP ES ASAT AS CLS PERSE RITA DTT CENT ADTRST ICI SEAT ERTC SCI RSTRS CRS CATIA RDIAIC LO TRICRT RET BOT ATT ATCO RI CRTC MIA ESATA TCR IT RECASCLS CROAT CIEL AT RET ATEASTAS TRUS TICS AICI TICECRSERS 
CHECKING RESULTS OF SECOND TEST SERIES ABOUT ACCEPTANCE 
IF (A.GT.3) THEN 
FLAG(1)=.FALSE. 
END IF 
IF (B.GT.3) THEN 
FLAG(2)=.FALSE. 
END IF 
IF (C.GT.5) THEN 
FLAG(3)=.FALSE. 
END IF 
IF (D.GT.3) THEN 
FLAG(4)=.FALSE. 
END IF 
IF(.NOT. (FLAG(1).AND.FLAG(2).AND.FLAG(3).AND.FLAG(4))) THEN 
FLAG(1)=.TRUE. 
FLAG(2)=. TRUE. 
FLAG(3)=.TRUE. 
FLAG(4)=. TRUE. 
GOnTO™ 142 
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Es 


+ 


+ 


110 


END IF 


alealeats leeks elo alraleale alates outetooies 
KKRKRKEREKRERERAERRERRREEREAEERREE BREE RORAAA SREB RREREREEREEREEEEEEEERER 


AFTER SECOND INSPECTION DETERMINATION OF LCB 
LCB=SELCB(A,B,C,D) 
HRA EREKEK ERR ERAERRERL RRR EERE R ERE REE ERERR EERE REREREEREEREE 
CHECKING FOR BONUS AFTER SECOND INSPECTION 
IF ( LCB.GT.BLINE ) THEN 
BSI=BSI+1 
GO TO 112 
END IF 
KERKEKRAEKKRAIRERRER ERE REE RR EE RE EERE RRR TERRE EERE ERR RRERREREER 
COUNTING FOR CHECKING 1000 REPLICATIONS 
COUNT = COUNT + 1.0 
A=0.0 
B=0.0 
C=0.0 
D=0.0 


ARAL EL IVIL IL EL AD FEIT ID GEIL AE GY AE AU AE IN IN EU IEEE ALIN GN ID IE IE ILD AD IU IN OR IR AN ER OE IS AVIV IVIL GL AIV GL GD IL ATID IL ER GE GL IL GL ad Gee ARR IN ED 


CHECKING FOR 1000 REPLICATIONS 

IF (COUNT.GT.1000) THEN 

GO TO 130 

END IF 
CONTINUE 
REPRE AIAREREERITEE EEE ERERERERER EIR ERATE EREECEREERERERERERERR ERE 
PERCENTAGE ESTIMATION OF GETING BONUS FOR FIRM A 
BTOT = BFI+BSI 
PRCT=BTOT/1000.0 


RKREKKEI SREERERERERIEERRIERERERERERERRERREREERERERCER TREE 


WRITING RESULTS 
WRITE (15,1) 
FORMAT (//,16X,' BONUS PLAN SIMULATION FOR FIRM A ',2X) 
WRITE, (15,2) 


FORMAT (16X, tttinnninidiinninkiniininciiiik | 2X, //) 


WRITE (15,3) 
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12 


13 


14 


FORMAT(4X,' FIRM A HAS FOLLOWING LONG RUN PROBABILITIES IN TESTS ' 
+,2X) 

WRITE (15,4) 

FORMAT ( 4X : TERA LEEK AES ARER RE ER ERE R REAR EER EEE BER EE EERE RE EREREERRERS ! 
tis 2X) 

WRITE (15,5) PSIM,PSITH,PSIV,PSIA 

FORMAT(4X, ‘PROBABILITY OF SUCCESS IN MANUFACTURER TEST IS',2X, 
+F8.6,//,4X,' PROBABILITY OF SUCCESS IN TEMP. AND HUM. TEST IS',2x, 
+F8.6,//,4X, ‘PROBABILITY OF SUCCESS IN VIBRATION IS',2x, 
+F8.6,//,4X, ‘PROBABILITY OF SUCCESS IN ALTITUDE IS',F8.6,2X,//) 
WRITE (15,6) BFI 

FORMAT(4X,'FIRM A GAT BONUS AFTER FIRST INSPECTION ',2X,F6.1,2X, 
+'TIMES' ,2X) 

WRITE (15,7) 

FORMAT ( 4X ; VRE ARERR RRA EERE RRA AR ARE ERS ARE RARER BREE ARERR ERRERES 
FRREEEK' IX, /) 

WRITE (15,8) BSI 

FORMAT(4X,'FIRM A GAT BONUS AFTER SECOND INSPECTION ',2X,F6.1,2X, 
+'TIMES' ,2X) 

WRITE (15,9) 


FORMAT (4X , | RRR RERARAR AAA ARERR ER ER ERE R ER EERE REE RERER ER ERRE 
Rake 2X, /) 

WRITE (15,11) BTOT 

FORMAT(4X, ‘TOTALLY FIRM A GAT BONUS IN 1000 REPLICATIONS ', 2x, 
oF5.1,2X, TIMES’ , 2X) 


WRITE (15,12) 
FORMAT ( 4X , * drtttirtteintedeat arise deat tetededted teeter k ERIKA RII ERR 


fre! 2X, /) 

WRITE (15,13) PRCT 

FORMAT(4X, ‘GETTING BONUS PERCENTAGE OF FIRM A IS',1X,F6.3,2X) 
WRITE (C15 14) 

FORMAT (4X, * ere A AAA ARR EER ARERR ERR ERE RRRERE RRR ER ARE 
pee ek) 


WRITE (15,15) BLINE 


13) 


15 FORMAT(4X,' BONUS LINE FOR FIRMS IS',1X,F6.4,2X) 
WRITE (15, 16) 
16 FORMAT (4X, | eeitreerarcce cree kk AAA RARHAA AAA RRRAARRK RARER AER ER RRR E RE 
+e OX ///) 
STOP 
END 
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APPENDIX I. 95 % LCB’S FOR DSBS ( EQUAL PROBABILITIES ) 


Table 18. 95 % LCB’S FOR DOUBLE SAMPLING BONUS SYSTEM 


95 % LCB See 95 % LCB 
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Table 19. 95 % LCB’S FOR DOUBLE SAMPLING BONUS SYSTEM 
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Table 20. 95 % LCB’S FOR DOUBLE SAMPLING BONUS SYSTEM 


FAILURE F FAILURE Sp 
ae 95 % LCB VECTOR 95 % LCB 


ann 
aia) | aones [i] 
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APPENDIX J. DOUBLE SAMPLING BONUS SYSTEM WITH EQUAL 


PROBABILITIES 


Table 21. DSBS (EQUAL PROBABILITIES) LCBFB = 0.825 


LOWER CONFIDENCE BOUND FOR BONUS IS 0.825 
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EQUAL TEST PROBABILITIES; LCBFB = 0.825 
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0.96 
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Figure 5. Double Sampling Bonus System With LCBFB = 0.825 
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Bonus percentages are tabulated and plotted below with LCBFB 0.850 


Table 22. DSBS (EQUAL PROBABILITIES) LCBFB = 0.850 


LOWER CONFIDENCE BOUND FOR BONUS IS 0.850 
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DOUBLE SAMPLING BONUS SYSTEM 
EQUAL TEST PROBABILITIES; LCBFB = 0.850 
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Figure 6. Double Sampling Bonus System With LCBFB = 0.850 
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Bonus percentages are tabulated and plotted below with LCBFB 0.875 


Table 23. DSBS (EQUAL PROBABILITIES) LCBFB = 0.875 


LOWER CONFIDENCE BOUND FOR BONUS IS 0.875 
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DOUBLE SAMPLING BONUS SYSTEM 
EQUAL TEST PROBABILITIES; LCBFB = 0.875 
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Figure 7. Double Sampling Bonus System With LCBFB = 0.875 
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Bonus percentages are tabulated and plotted below with LCBFB 0.900 


Table 24. DSBS (EQUAL PROBABILITIES) LCBFB = 0.900 


LOWER CONFIDENCE BOUND FOR BONUS IS 0.900 
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DOUBLE SAMPLING BONUS SYSTEM 
EQUAL TEST PROBABILMES; LCBFB = 0.900 
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Figure 8. Double Sampling Bonus System With LCBFB = 0.900 
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Bonus percentages are tabulated and ploted below with LCBFB 0.950 


Table 25. DSBS (EQUAL PROBABILITIES) LCBFB = 0.950 


LOWER CONFIDENCE BOUND FOR BONUS IS 0.950 
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DOUBLE SAMPLING BONUS SYSTEM 
EQUAL TEST PROBABILITIES; LCBFB = 0.950 
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Figure 9. Double Sampling Bonus System With LCBFB = 0.950 
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Bonus percentages are tabulated and ploted below with LCBFB 0.999 


Table 26. DSBS (EQUAL PROBABILITIES) LCBFB = 0.999 


LOWER CONFIDENCE BOUND FOR BONUS IS 0.999 
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DOUBLE SAMPLING BONUS SYSTEM 
EQUAL TEST PROBABILMES; LCBFB = 0.999 
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Figure 10. Double Sampling Bonus System With LCBFB = 0.999 
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APPENDIX kK. 95 % LCB’S FOR DSBS (DIFFERENT PROBABILITIES) 
Bonus percentages are tabulated and ploted with different probabilities. 
LCBFB:= 0825 


Table 27. DSBS (DIFFERENT PROBABILITIES) LCBFB = 0.825 


DOUBLE SAMPLING BONUS SYSTEM 
LCBFB = 0.825 


Figure 11. Double Sampling Bonus System With LCBFB = 0.825 
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Bonus percentages are tabulated and ploted with different probabilities. 
LCBFB = 0.850 


Table 28. DSBS (DIFFERENT PROBABILITIES) LCBFB = 0.850 


DOUBLE SAMPLING BONUS SYSTEM 
LCBFB = 0.850 
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Figure 12. Double Sampling Bonus System With LCBFB = 0.850 


143 


Bonus percentages are tabulated and ploted with different probabilities. 
LCBFB = Ui57/5 


Table 29. DSBS (DIFFERENT PROBABILITIES) LCBFB = 0.875 


DOUBLE SAMPLING BONUS SYSTEM 
LCBFB = 0.875 
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Figure 13. Double Sampling Bonus System With LCBFB = 0.875 


Bonus percentages are tabulated and ploted with different probabilities. 
LCBFB = 0.900 


Table 30. DSBS (DIFFERENT PROBABILITIES) LCBFB = 0.900 


DOUBLE SAMPLING BONUS SYSTEM 
LCBFB = 0.900 
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Figure 14. Double Sampling Bonus System With LCBFB = 0.900 
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Bonus percentages are tabulated and ploted with different probabilities. 
LCBFB = 0.950 


Table 31. DSBS (DIFFERENT PROBABILITIES) LCBFB = 0.950 


DOUBLE SAMPLING BONUS SYSTEM 
LCBFB = 0.950 
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Figure 15. Double Sampling Bonus System With LCBFB = 0.950 
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Bonus percentages are tabulated and ploted with different probabilities. 
LCBFB = 0.999 


Table 32. DSBS (DIFFERENT PROBABILITIES) LCBFB = 0.999 


DOUBLE SAMPLING BONUS SYSTEM 
LCBFB = 0.999 
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Figure 16. Double Sampling Bonus System With LCBFB = 0.999 
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